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ABSTRACT 


The transformation of surface gravity waves propagating through shallow regions 
is investigated with extensive field data from the North Carolina continental shelf. A spec¬ 
tral energy balance equation is derived for a bidimensional bottom topography with ran¬ 
dom small-scale irregularities, in which bottom friction is introduced heuristically whith 
a parameterized source term, and solved numerically using a hybrid Eulerian-Lagrangian 
scheme. This new model named CREST (Coupled Rays with Eulerian Source Terms) deter¬ 
mines accurately refraction of waves by subgrid-scale depths variations using precomputed 
rays, allowing applications to large coastal areas with relatively coarse grids. Hindcasts 
of swell events during field experiments show large variations in wave heights caused by 
refraction and bottom friction. Widespread observations of sand ripples confirm that the 
bottom roughness is enhanced by wave-generated vortex ripples, thus sheltering the shore 
from offshore swells by dissipating wave energy in the bottom boundary layer. Resulting 
wave height attenuation up to 70 % (84 % of the wave energy) was observed in moder¬ 
ately energetic conditions. Bragg scattering of waves by wavelength-scale bottom features 
significantly increases (up to a factor two) the directional spread of waves. 
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RESUME 

La transformation d’ ondes de gravite de surface est etudiee au moyen d’ une equa¬ 
tion d’ evolution pour les densites spectro-angulaires, etablie a partir des equations de con¬ 
servation de la quantite de mouvement et de la masse, dans le cas d’ un fond lisse avec une 
topographie irreguliere. L’effet de la rugosite du fond est introduit de maniere empirique. 
Un modele numerique hybride Eulerien-Lagrangien est presente pour resoudre cette equa¬ 
tion. Entre autres avantages ce nouveau modele, denome Cretes (Couplage de Rayons par 
Termes de Source Euleriens) permet d’utiliser des maillages non-structures pour le terme 
de source, tandis que les effets de la refraction et du levage causes par les variations de 
profondeur aux echelles sous-maille sont precisement pris en compte grace aux rayons pre- 
calcules. Ainsi ce modele pent etre applique a de vastes zones cotieres avec des maillages 
relativement laches. 

Des observations de rides de sable sur le plateau continental de Caroline du Nord, 
combinees a des simulations numeriques d’ episodes houleux observes au cours de cam- 
pagnes a la mer, montrent que les variations de la hauteur des vagues a la traversee d’ un 
plateau continental large, peu profond et sableux, sont essentiellement causes par la refrac¬ 
tion et la friction au fond, pour laquelle les rides crees par les vagues jouent un role deter¬ 
minant. Ces rides protegent la cote des houles d’ amplitude moyenne, avec, pour les cas 
decrits ici, une attenuation maximale de 70 % de la hauteur de vague, soit 84 % de I’energie. 

Ees proprietes directionnelles de la houle sont modifiees par la presence d’ irregu- 
larites bathymetriques d’ echelle horizontale comparable a la longueur d’ onde de la houle. 
Cette interaction houle-topographie est decrite ici par la theorie de la diffusion de Bragg 
pour des vagues aleatoires, et representee par un terme de source, ou plus precisement 
d’ echange, dans I’equation d’ evolution des densites spectro-angulaires. Sur le plateau 
Nord-Carolinien cette diffusion pent augmenter I’etalement angulaire des vagues d’ un 
facteur deux, contrariant le retrecissement du spectre angulaire cause par la refraction. Ces 
previsions correspondent aux observations, malgre la persistance d’ une erreur residuelle 
qui pent etre attribute en grande partie aux incertitudes sur les conditions a la limite au 
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large, mais qui peut aussi etre la trace d’ autres phenomenes non prepresentes dans le 
modele. 
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I. INTRODUCTION 

A. MOTIVATION AND OBJECTIVES 

Surface gravity waves, in oceans, lakes, or ponds, often simply called ‘waves’, are 
one of the earliest observed and studied among the many phenomena that physicists refer 
to as waves. Waves inspired many of the greatest mathematicians, physicists, and oceanog¬ 
raphers from the 19* and 20* centuries, from Boussinesq, Laplace, Stokes, and Rayleigh, 
to Hasselmann, Munk, Phillips, Sverdrup, and Zakharov, to name only a few. Because of 
their simple yet rich nature (waves are generally dispersive and weakly non-linear), and 
the natural homogeneity of water bodies and the atmosphere above, waves were the natural 
testing ground for many theories and ideas, from Fourier analysis, initially developed for 
heat conduction, to Fermi-Pasta-Ulam recurrences or weak turbulence theory. In the view 
of this enormous body of knowledge, one may think that everything about waves is already 
perfectly known, but this is not the case. In particular the generation of waves by wind, and 
the resulting highly non-linear breaking waves are still the subject of active research. As 
waves propagate away from the area where they were generated by wind, they evolve into 
regular swell that can radiate across large ocean basins with little energy loss, until they 
reach the coast and finally dissipate in the surf zone, or reflect from a cliff. Both seabed ge¬ 
ology and topography strongly influence the evolution of swell, across continental shelves, 
shallow marginal seas, or lakes, in a way that is still poorly understood. 

The goal of the present dissertation is to determine these effects quantitatively, from 
theory and wave measurements in the field, using a numerical model. The mathematical 
model and adapted numerical techniques to represent swell evolution proposed here will 
hopefully not only improve conceptual understanding of the physical processes, but also 
lead to more accurate predictions and hindcasts of swells near the shore. 

In the marine environment all human activities depend on such wave predictions, 
from fishing, shipping and naval or amphibious operations, to the development of coastal 
regions and the offshore industry. Preparation for D-day during the second world war was 
the driving force behind some of the first attempts at systematic wave prediction using 
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available theory (Sverdrup & Munk, 1947). Just like waves are driven by the wind, so were 
wave models pushed by the advent of operational weather forecasting and its refinements 
over the last sixty years. Observations and predictions of rising sea level and storm fre¬ 
quency, both associated with global warming, are increasing the challenge posed by ocean 
waves, and continue to call for better wave predictions. Current wave models used for engi¬ 
neering design or operational predictions already provide immensely valuable information, 
but some important physical processes at play in the evolution of waves in shallow water 
are still poorly known. Models therefore rely on empirical relations that are derived from 
limited data sets, and thus cannot be use with confidence in untested situations. 

The dissipation of wave energy due to friction on a sandy bottom offers an inter¬ 
esting example. Based on the ideas used for ocean currents and air flows, Hasselmann 
& Collins (1968) proposed that this bottom friction may follow a quadratic drag law and 
they derived the corresponding effect for random waves. Data from the 1968 JOint North 
Sea Project (JONSWAP, Hasselmann et al, 1973), off the coast of Sylt in the North Sea, 
provided a real-life test of their theory. The theory failed to explain the observed wave 
dissipation variations, even after accounting for the dominant tidal currents. Lacking a 
reliable physics-based bottom drag formulation, most wave models today use a simple pa¬ 
rameterization of bottom friction based on the quadratic drag law for waves in the presence 
of dominant tidal currents, with a tunable coefficient. This coefficient is usually set at the 
average measured value during JONSWAP, although its variations span two decades in the 
JONSWAP dataset. It is a tribute to human engineering genius that this parameterization 
generally gives reasonable results, even in the absence of significant currents. However 
this crude representation of bottom friction certainly lacks the reliability and accuracy of 
theories based on first principles foundation on which one may like to build a multi-billion 
dollar harbor, or plan for the next D-Day. Fortunately a better understanding of the bottom 
friction is now available, as the thin bottom boundary layer where energy dissipation takes 
place has been studied theoretically, in the laboratory, in numerical experiments, and, to a 
lesser degree, in the field. However even the most elaborate boundary layer models (e.g. 
Reichardt, 1951; Weber \99\b) need some information on the roughness that depends on 
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the bottom geology and bedforms. The strong effect of wave generated sand ripples on 
bottom drag was already measured in 1963 in the laboratory (Zhukovets, 1963), but the 
introduction of this effect in wave forecasting models (Graber & Madsen 1988; Graber, 
Beardsley & Grant, 1989; Tolman, 1994) had to await the first parameterizations of mov¬ 
able sand bottom roughness by Grant & Madsen (1982), and Madsen et al. (1990). Before 
the present work no field data was available to demonstrate the importance of sand rip¬ 
ples for wave evolution and test these parameterizations. This slow progress should not 
overshadow successes in other areas of wave research, such as depth-induced breaking, 
but it illustrates how science may be hindered by the lack of good public data and the 
fragmentation of the research community into very small groups, despite international fed¬ 
erating efforts such as the meetings of the WAve Model Development and Implementation 
(WAMDI) group, and the Waves In Shallow Environments (WISE) group. 

The SHOaling Waves Experiment (SHOWEX), conducted in 1999 and funded by 
the U.S. Office of Naval Research, provided a much needed impetus to take a fresh look at 
the physics and the parameterizations of waves in shallow water. It provided a rich dataset 
of wave measurements, bathymetric and geological surveys, that is used in the present 
dissertation. The numerical wave model developed in chapter II was first validated with 
data from a previous experiment (DUCK94), showing the effect of sand ripples on bottom 
friction, but the theoretical results on wave scattering by irregular small scale topography 
(chapter III), could not have been tested without the high-resolution bathymetric surveys 
performed during SHOWEX. The presence of sand ripples, analyzed in chapter IV, was 
determined for the first time on a large scale by sidescan sonar surveys conducted during 
SHOWEX, providing a unique combination of wave measurements and bedform images. 
The wide range of wave conditions during this three-month long experiment, including 
major hurricanes and northeaster storms, allowed for a comprehensive statistical validation 
of the representation of sand ripples and wave-topography scattering in the wave model, 
in chapter V. Perspectives arising from that data set and the new numerical model are 
presented in chapter VI, with a summary of the main results of this dissertation. 
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Before getting into these new developments the non-expert reader may consult, for 
an in-depth introduction, the following overview of related wave theories, models, and 
observations. 

B. STATISTICAL DESCRIPTION OF WAVES 

Waves at the air-water interface modify the surface elevation, and the pressure and 
velocity below and above the surface. Wave motion is characterized by oscillations around 
a position of equilibrium, a ‘horizontal’ sea surface, under the action of gravity that tends 
to restore the equilibrium. In the presence of wind or steep wave fronts, surface tension 
may act as another restoring force for very short (capillary) waves (figure 56). This latter 
force will be neglected in the present work where the focus is on long period waves (fre¬ 
quency / less than 0.25 Hz), in the absence of wind. The position of the surface is defined 
by the function z = ^ (x, t) where x is the horizontal position vector, z is the vertical coor¬ 
dinate, pointing upward, with z = 0 the equilibrium position of the surface, and t is time. 
Assuming irrotational wave motion, a very good approximation in general, the horizontal 
velocity vector u and vertical velocity w can be expressed as gradients of a velocity poten¬ 
tial (|)(x,z,t): u = V(|), where V is the horizontal gradient operator, and w = d(|)/dz. The 
pressure perturbation p (x, z, t) can be given in terms of (|) and ^ using Bernoulli’s equation 
(e.g. § III.A). The wave motion is therefore completely described by (|) (x, z, t), which obeys 
Laplace’s equation, and ^(x,t). 

Since the objective of the present work is to describe the evolution of such waves 
over distances that are generally much larger than the wavelength, it is not practical to 
determine accurately ^ and (|) everywhere. Instead a statistical representation of the sea 
surface will be preferred (see Komen etal., 1994). The probability distribution of the sea 
surface elevation ^ at any horizontal location is generally very nearly Gaussian as the result 
of the superposition of many different waves with different frequencies and directions, 
propagating at different speeds, except in very shallow water where the wave height is of 
the same order as the water depth (Elgar & Guza, 1985). As a result of this property the 
joint probability distribution of the surface elevation can be determined from the surface 
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elevation covariance function 

(^(Xi,ti)^(x2,t2))- (I-l) 

A similar covariance can be used to describe the velocity potential. 

The three-dimensional Fourier transform of (LI) with respect to the two horizon¬ 
tal space and the time coordinates yields the (vector) wavenumber-frequency spectrum 
E (k,/). Since wave properties usually evolve slowly in space and time, the Fourier trans¬ 
form, a decomposition in sinusoidal waves, is advantageously replaced by a decomposition 
of ^ and (|) in slowly modulated sinusoidal waves, that may be based on the evolutionary 
spectral theory of Priestley (1965), giving a spectrum £ (k,/,x,t) that is also a slow func¬ 
tion of horizontal space and time, (see Komen et al, 1994, and chapter III). This spectrum 
can be evaluated by applying a three-dimensional Fourier Transform to time series of sur¬ 
face elevation maps, ^(x,t), or related quantities such as the surface slope, obtained by 
continuously imaging the sea surface in a small region with a dwelling aircraft (Dugan 
et al, 1996). 

A projection of this three-dimensional spectrum into two, or just one, dimensions is 
more readily accessible from simpler measurements. In linear theory, a very good approx¬ 
imation for waves in general and swell in particular, the wave frequency / and wavenum¬ 
ber magnitude A: = |k| are related (in the absence of currents) by the dispersion relation 
// (27t) = gA:tanh(A:/z), where g is the apparent gravity acceleration and h is the water 
depth. Therefore the projection of the three-dimensional spectrum into a two-dimensional 
wavenumber E (k) or frequency-direction E (/, 0) spectrum loses only the information on 
the weak non-linear effects. The wavenumber spectrum can be estimated, more or less 
directly and for different ranges of wavelengths depending on the platform, by surface 
imaging instruments such as scanning altimeters (e.g. Hwang et al, 2000), Real Aperture 
Radars (Jackson, Walton & Baker, 1985), and Synthetic Aperture Radars (e.g. Komen 
et al, 1994). In-situ measurements of velocity, pressure, or surface elevation with coherent 
arrays can provide estimates of E (/, 0) (e.g. Davis & Regier, 1977; Long & Hasselmann, 
1979; Pawka, 1983; Herbers & Guza, 1990). E (/, 0) is convenient for practical use as / is 
conserved during wave propagation (in the absence of currents and neglecting non-linear 
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effects) whereas k is modified by depth shoaling. An example spectrum is given in figure 1, 
using data from a coherent array of pressure sensors in 8 m depth at the U. S. Army Corps 
Field Research Facility in Duck, North Carolina. These estimates are more precise as one 
increases the number of sensors in the array. Most analysis methods for this type of data 
involve the cross and auto-correlations between the various measurements, assuming uni¬ 
form statistics over the length of individual records. Alternative analysis techniques have 
been proposed to give higher resolution in the directional spectra by reducing the frequency 
resolution. Such a technique developed by Donelan, Drennan & Magnusson (1996) gives 
interesting results but interpretation is made difficult by the unfounded hypothesis that for 
any given time and frequency band the wave field is dominated by one wave group coming 
from a well-defined direction. 

The projection of E (/, 0) on the frequency axis (the frequency spectrum, see figure 
lb) can be evaluated from the time series of surface elevation (as measured by an infrared 
laser mounted on a tower) or pressure at a fixed sub-surface elevation, for example at the 
sea floor. It is the most common in-situ measure of the sea state but it gives no information 
on the directional properties of waves. 

The main wave properties are often quantified with a few parameters. The most 
common are the significant wave height Hs, defined as four times the square root of the 
elevation variance (the area under the curve in figure l.b), the peak frequency fp at which 
the spectrum E (/) is maximum, or its reciprocal, the peak period Tp = 1/ fp. In the limit 
of a narrow frequency spectrum, Hs is equal to the average height Hi /3 of the highest one 
third of the waves. Hs is indirectly determined from the shape of satellite radar altimeter 
pulses (e.g. Rao et al, 1990), back-scattered by the sea surface, providing the only wave 
measurement available globally that is used for operational wave forecasting. 

The elevation variance density spectra E defined above in wavenumber-frequency, 
wavenumber, frequency-direction, or frequency spaces are commonly called wave energy 
spectra. However they must be multiplied by the water density p and gravity g to have units 
of energy. For the example spectrum shown in figure l.a, the onshore energy flux, that is 
the average rate (over many waves) at which wave energy is released on the nearby beach 


6 





0 19 38 57 76 95 

Energy density (in- s^) 



frequency (Hz) 

Figiue 1. Wave spectra 

(a) Example fr equency-direction wave spectnim, determined fr om pressure aixay measure¬ 
ments in 8 m depth at Duck, NC, October 19, 1994, 7:00-10:00 EST. (b) Coiiesponding 
frequency spectnmi. in which the first hannonic peak at f = If pis probably due to non¬ 
linear effects, which are enhanced in shallow water and for large amplitude waves such as 
these (see § lI.E. 1 for more infoiination on this event). 
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per meter of coastline, is 13 kW. This estimate illustrates how powerful waves are, even 
after losing 75% of their energy across the shelf because of bottom friction (see chapters II 
and V). 

Between coherent arrays providing detailed directional information, and single sen¬ 
sors giving simple frequency spectra, intermediate self-contained devices have been devel¬ 
oped to measure some representative directional properties of waves. The directional prop¬ 
erties of the wave spectrum are usually summarized with a mean direction for the entire 
spectrum, 0, for each frequency, 0 (/), or at the peak only, 0^, as indicated on figure 1. The 
corresponding directional spread Oe, o© (/), or Oe^p, gives a measure of the half-width of 
the directional distribution of the wave energy (Kuik, van Vledder & Holthuijsen, 1988). 
In the present dissertation Oe is equal to the standard deviation in radians if the directional 
distribution is narrow, and reaches a maximum value of 2^/^ radians (that is 81 degrees) 
for an isotropic distribution. These simple parameters can be determined from the first 
Fourier components of the directional distribution at each frequency. Both the first and sec¬ 
ond Fourier components can be estimated from collocated tri-axis acceleration time series 
(as measured by a Directional Waverider buoy), or any combination of three independent 
scalar quantities related to these by linear wave theory, such as the heave, pitch, and roll of 
a floating platform (Longuet-Higgins, Cartwright & Smith, 1963; Long, 1980), or pressure 
p and horizontal velocity components u and v. A ‘cloverleaf’ curvature buoy (Cartwright 
& Smith, 1964), was designed to provide more Fourier components, but proved delicate to 
use. 

Moored surface-following buoys can be used in any water depth for accurate mea¬ 
surement of the dominant swells (e.g. O’Reilly et al, 1996), but these large Lagrangian 
sensors do not resolve short, nonlinear waves, and are suspected to miss extremely large 
(so-called ‘freak’) waves by sliding off their top. Bottom-mounted pressure and p-u-v 
gauges are more accurate but limited to shallow areas (~ 20 m, depending on the fre¬ 
quency of interest) because of the vertical exponential attenuation of the wave signal over 
a height equal to the wavelength. 

The directional information given by the first two Fourier components measured 
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by all these instruments can be used to estimate a full frequency-direction spectrum. The 
Maximum Entropy Method (Lygre & Krogstad, 1986), used in the present work, gives an 
estimate that is constrained to fit the data exactly but the construction of a directional dis¬ 
tribution from only four parameters calls for caution in interpreting the resulting estimates. 


C. THE ENERGY BALANCE EQUATION 

The evolution of the wave spectrum E (k,x,t), can be described with a spectral 
energy balance (Gelci, Cazale & Vassal, 1957). Neglecting the effects of currents, it is 
given by (e.g. Whitham 1974): 

dE 

-^ + ■ (CgE) + Vk ■ (CkE) = S (1.2) 


where Vx and Vk are horizontal divergence operators in geographical and wavenumber 
space respectively, and Cg (the group speed) and are the corresponding energy transport 
velocities. The source term 5(k,x,t) is the net rate of energy transfer to component k 
resulting from interactions between wave components, or interactions with the bottom, 
the atmosphere, or other flows in the water column such as internal waves. The effects of 
currents can be added by replacing the energy spectrum by the action spectrum N = 2nE //, 
and replacing S with the corresponding action source term (e.g. Bretherton and Garret, 
1969; Komen et al, 1994). 

The spectral energy balance can also be formulated from a Lagrangian point of 


view: 



(1.3) 


where the left-hand side is the rate of change of E following a wave component along its 
ray trajectory. (1.2) and (1.3) are equivalent for waves that obey a dispersion relation of the 
form / = Vk (k, x, t) where / is the wave component frequency. However, in contrast to 
the Eulerian balance (1.2), the along-ray conservation of spectral densities for 5 = 0 is valid 
only in k-space (Eonguet-Higgins, 1957). In both cases the left hand sides represent the ad- 
vection of wave energy accounting for refraction (the turning of wave crests around shoals 
and depressions) and shoaling (the increase in energy when the group speed decreases). 

9 



The stochastic representation outlined above, generally referred to as weakly non¬ 
linear and ‘phase averaged’, assumes that phases of spectral wave components are uncor¬ 
related. This is appropriate for waves in intermediate and deep water where nonlinearity is 
weak and phases are randomized by dispersion (Benney & Saffman, 1966). This simplifica¬ 
tion is not applicable in very shallow water, inside or very near the surf zone, where waves 
are weakly dispersive and non-linearity is enhanced. The reader may follow Freilich & 
Guza (1984) and Herbers & Burton (1997) for spectral approximations used in this regime, 
or Osborne et al. (1998), for a different approach based on inverse-scattering theory. 

The energy balance (1.2) or (1.3) provides a simple prognostic equation for the evo¬ 
lution of the wave spectrum. In shallow to intermediate water depths, surface gravity waves 
are affected by sea bed features with a wide range of scales (figure 2). Wave refraction over 
large scale (nominally 1 to 10 km) bottom features can induce dramatic variations in wave 
energy along the coast that are readily observed (e.g. Munk and Traylor, 1947). The effects 
of refraction on the evolution of wave spectra are generally well understood, and accurately 
represented in the left-hand sides of (I.2-I.3) by geometrical optics models (e.g. Longuet- 
Higgins, 1957; Mei, 1989; O’Reilly & Guza, 1993; §II.A). 

However the ray trajectories followed by wave groups, predicted by the geometrical 
optics approximation, may be modified by the strong non-linearity of steep waves (Wille- 
brand, 1975). This effect will be neglected in the present work as it is weak for the low 
steepness swell observed outside storm wave-generation regions. The present analysis is 
focused on swell, with small surface slopes in light wind conditions, so that generation of 
waves by the wind, dissipation of energy due to wave breaking, and mutual interaction of 
different wave trains can be neglected. These processes that control wave evolution in deep 
water are usually represented in the right-hand side source terms of (I.2-I.3), by the sum of 
three parameterized source terms. Sin for the wind energy input, Sdis for the wave breaking 
dissipation, and Sni for the non-linear wave-wave interactions. Although the latter is never 
completely absent even without wind, it can be discarded here as it modifies the wave spec¬ 
trum on a time scale of the order of (e^/p) ^ where £ is the mean wave slope, and fp is 
the peak frequency (Hasselmann, 1962). This time scale for the lowest order interactions 
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Figiue 2. Wave-bottom iuteractions 

Smnmaiy of the effects of different bottom imdulation scales (the bottom .r-axis coordi¬ 
nate 27t// is the reciprocal bottom wavelength), for typical moderate swell motion scales 
(wavelength kjln and significant near-bed orbital displacement diameter d^, indicated by 
vertical dash-dotted lines). The thick curve is a typical bottom slope spectnmi for the North 
Carolina shelf derived fiom bathymetry siuveys for bottom wavelengths larger than 40 rn. 
A variety of bedforins of wavelengths less than 10 m can be generated by the waves de¬ 
pending on the forcing conditions. Of particular interest are steep and active vortex ripples 
(sohd cmwe, with wavelengths and heiglits estimated fiom sidescan sonar siuveys) main¬ 
tained by the wave orbital motion. They str ongly enhance the tmbuleut dissipation of wave 
energy on the shelf These ripples become relic in benign wave conditions, or are oblit¬ 
erated in strong forcing conditions, timriug the sea bed into a sheet flow layer with lower 
drag (dashed ciuve). 
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is generally much longer (of the order of 20 days for 0.1 Hz swell with £ = 0.05) than the 
propagation time of swell over the regions of interest in the present work (about three hours 
for 0.1 Hz swell crossing the North Carolina continental shelf). Higher order interactions 
are even weaker (Krasitskii, 1994). 

Swell evolution is therefore determined primarily by wave-bottom interaction pro¬ 
cesses. Small scale (one half to several wavelengths) bottom features can scatter waves, 
and bottom irregularities shorter than the typical orbital diameter of the wave motion at the 
bottom, can be considered as roughness elements that are responsible for the dissipation of 
wave energy in the bottom boundary layer (figure 2). These latter two effects are discussed 
in the next two sections together with their representation as source terms in the right hand 
sides of (1.2-1.3). 

D. WAVE-BOTTOM SCATTERING 

Hasselmann (1966) proposed a statistical theory for the evolution of random sur¬ 
face gravity waves over an irregular bottom assuming spatially homogeneous conditions 
(i.e. uniform surface wave and bottom elevation spectra). At the lowest order, two wave 
components with the same radian frequency (0 but different wavenumber vectors k and k' 
exchange energy in a resonant triad interaction with the bottom component that has the 
difference wavenumber 1 = k — k' (figure 3). This Bragg scattering process is potentially 
important for the directional properties of the waves. Long (1973) applied Hasselmann’s 
theory to swell in the North Sea with some assumptions about the unknown statistical prop¬ 
erties of the bottom topography. His results suggested that back scattering of surface waves 
from bottom undulations with wavelengths close to half the surface wavelength (k —k', 
1 ^ 2k) could explain the swell energy decay observed during the JONSWAP experiment 
(Hasselmann etal., 1973), but subsequent bathymetric surveys (Richter, Schmalfeldt & 
Siebert, 1976) showed that the amplitude of seabed undulations at the site of the JON¬ 
SWAP experiment was too small to cause significant back scattering. Ewing & Pitt (1982) 
observed reflected waves increasing in height away from a rocky coast, qualitatively consis¬ 
tent with wave-bottom Bragg scattering. Yet the importance of wave scattering by natural 
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Figiue 3. Resonance condition for Bragg scattering 
The interaction between a siuface wave with wavenumber k' and a bottom component with 
wavemunber 1 excites a smface wave with the siun wavenumber k = k' + 1 . For fixed k , 
the resonant k' and 1 lie on the solid and dashed circles, respectively. 

seabeds has remained imknown in general. It is believed to be weak, and the lack of detailed 
bathymetric data has prevented fiulher investigations (Komen et ah, 1994). 

Using a different detemirnistic approach, Davies (1979) derived an analytical solu¬ 
tion for the weak reflection of a monochr omatic wave train propagating at normal incidence 
over a patch of simrsoidal bars, that was srrbseqirently verified in laboratory experiments 
(Heathershaw 1982; Davies & Heathershaw, 1984). Davies’ theory does not accoimt for 
the decay of the incident wave, losing energy to the reflected component, and therefore 
overestimates strong reflections, in particirlar at resonance where 1 = 2k. Mei (1985) de¬ 
rived a more acciuate energy conserving sohrtion, valid close to resonance, that was con¬ 
firmed by experiments (Kara & Mei, 1987). The more general case of obliqite incidence 
was considered by Mei (1985), Dahymple & Kirby (1986) and Kirby (1993). Mei (1985) 
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further generalized his theory to bars superimposed on a sloping bottom. Kirby (1986a; 
1986Z?) subsequently showed that Mei’s (1985) generalized theory ean also be derived from 
modified mild slope equations. Other related developments inelude non-linear effeets in a 
long wave approximation (Benjamin, Boezar-Karakiewiez & Pritehard, 1987), higher or¬ 
der Bragg seattering (Mitra & Greenberg, 1984; Belzons, Rey & Guazzelli, 1991; Liu & 
Yue, 1998; Agnon & Sheremet, 2000), extended mild slope equations for steep topography 
(Athanassoulis & Belibakis, 1999), and investigations of Anderson loealization of waves 
on a random bottom (Devillard, Dunlop & Souillard, 1988; Belzons, Guazzelli & Parodi, 
1988). Implieations for sediment transport and the formation of multiple sand bar systems 
just outside the surf zone were diseussed by Heathershaw (1982) and Mei (1985). 

Whereas Hasselmann’s (1966) stoehastie theory gives an energy balanee equation 
that is an effieient tool for predieting the speetral evolution of random waves, it is restrieted 
to homogeneous wave and bottom topography properties, and has not been verified ex¬ 
perimentally. In eontrast, Mei’s (1985) deterministie theory is more general and has been 
verified for simple eases, but it has not been applied yet to a natural sea bed beeause it 
requires a numerieal solution to an elliptie equation that is prohibitively expensive for large 
domains. Kirby (1986a) diseussed these two eomplementary theories but eould not reeon- 
eile them for the ease of monoehromatie waves traveling over a sinusoidal bottom. Indeed, 
Hasselmann’s theory assumes that the wave energy speetrum is eontinuous aeross the res- 
onanee manifold in order to determine its long-term evolution, and thus eannot be applied 
to monoehromatie waves (see Hasselmann, 1962, and Komen et al, 1994, for detailed dis- 
eussions of the eontinuum approximation in random wave seattering theory). 

In ehapter III we examine the effeets of wave-bottom seattering from natural seabed 
topography by extending Hasselmann’s (1966) theory to a heterogeneous wave field and 
sea bed topography. Hasselmann’s seattering souree term is re-derived on a gently sloping 
bottom with slowly varying wave and bottom speetral properties, eorreeting for an appar¬ 
ent error in the wave-bottom eoupling eoeffieients given by Hasselmann (1966) and Long 
(1973), and further redueing the likely importanee of Bragg seattering during JONSWAR 
In § III.B this result is verified through eomparisons with analytie results of deterministie 
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theories for waves propagating over a finite patch of sinusoidal bars, and the consequences 
of introducing a large-scale bottom slope are investigated. The effects of Bragg scattering 
on swell propagation across the North Carolina continental shelf are illustrated in § III.C 
with an implementation of the scattering source term in the spectral wave prediction model 
CREST (chapter II) using measured wave spectra and high-resolution bathymetric data. 
Extensive comparisons with field data are performed in chapter V. 


E. ENERGY DISSIPATION IN THE BOTTOM BOUNDARY LAYER 
1. Bottom drag 


In addition to the relatively well understood energy conserving wave-bottom inter¬ 
action processes, wave evolution across continental shelves and in shallow marginal seas is 
also believed to be strongly affected by non-conservative bottom boundary layer processes 
(e.g. Shemdin et al, 1980; Bouws and Komen, 1983; Weber, 1988; Young & Gorman, 
1995; Herbers, Hendrickson & O’Reilly, 2000). 

Eor sea beds composed of non-cohesive sandy sediments, the dissipation of wave 
energy in the bottom boundary layer was shown, in laboratory experiments, to be strongly 
dependent on the presence of sand ripples formed by the near-bed wave orbital motion (e.g. 
Zhukovets, 1963; Nielsen, 1992). Neglecting currents unrelated to the waves, the bottom 
boundary layer can be classified in three regimes, based on the ratio of friction and buoyant 
forces acting on a sand grain, and represented by the maximum Shields number \|/max (often 
denoted e^^ax) 


Vmax — 


f' 

J w“max 


(1.4) 


(5-l)gZ)’ 

where /(( is a skin friction factor, s is the specific density of the sand grains, D is their 
diameter, and g is the acceleration of gravity (Shields, 1936). Eor small values of V[t max , 
the bottom morphology does not change, thus retaining the history of past wave events and 
biological activity. In this ‘relic roughness’ regime wave energy dissipation is minimal as 
bottom velocities are small and turbulence is weak. 

In the following the magnitude of the wave energy dissipation will be measured by 
the dissipation factor fg, an averaged drag coefficient that gives the energy dissipation from 
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Figiue 4. The thiee boundary layer regimes 

(a) relic roughness, (b) active ripples, (c) sheet flow. In each the wavelength and water 
depth should be much larger than indicated, typically 100 m, versus 1 m for the bottom 
ripples. The horizontal aiTows indicate the velocity profile imder the wave crests (the flow 
reverses imder the troughs). Cmwed anows in (b) represent vortices generated in the lee of 
the ripple crests, to the right imder wave crests and to the left imder wave tiouglis. 


a representative bottom velocity cubed, /g mcreases dramatically as siuftcial sediments, 
intennittently set in motion by the wave orbital flow, organize into ripples that increase 
bottom roughness. This transition occms when v(/max increases past a threshold value V|/c 
(typically 0.03-0.1 for well-sorted quartz sand) (e.g. Nielsen, 1981). These ‘active ripples’ 
sharply increase the tiubulent dissipation of wave energy, as vortices are shed by the orbital 
flow at the ripple crests. According to experiments by Madsen, Mathisen & Rosengaus 
(1990) with random waves, /g is maximirm when V|/nns — l-2v|/c, where V|/niis is computed 
from the root mean sqirare bottom velocity amplitude usmg a formulation similar" to (1.4). 
Example swell conditions with this maximmn drag are a peak period Tp-\2 s, and signif¬ 
icant wave height Hs= 1.5 m in 25 m depth over well sorted quartz sand with grain size 
Z) = 0.15 mm. For larger values of V|/max (larger wave height or frequency), /g gradually 
decreases as ripple crests are eroded by stronger flow. 

For ver"y large values of i|/niaxn of the order of 10i|/g (Li and Amos, 1999), a layer 
of sediment, called ‘sheet flow’, moves with the water cohrmn, washing out ripples, giving 
a relatively weaker energy dissipation in the bottom boimdar'y layer. Both the thickness 
of this ‘sheet flow’ layer and the dissipation factor /g increase with i|/. These three flow 
regimes, relic roughness, active ripples, and sheet flow, are illustrated in figme 4. 


16 













































2. Bedform dynamics 

Since small-scale bedforms on the sandy sea floor play an important role in the 
transformation of waves aeross the shelf, we give here a brief review of studies on their 
properties. Bedforms not only affeet waves but also have a strong impaet on the transport of 
sediments, either as a result of their migration or beeause of their influence on the flow that 
shapes them. Their presence, formation, and evolution have been observed extensively in 
nearshore environments (e.g. Hunt, 1882; Forel, 1888; Dingier, 1974; Vineent & Osborne, 
1993; Gallagher, Elgar & Thornton, 1998; Traykovski et al, 1999). In the absence of mean 
eurrents, waves can generate ripples that are symmetric in cross-section. The formation of 
such wave ripples was first investigated in the laboratory by Darwin (1883) using a rotating 
bath. He noted the important role of the vortiees generated in the lee of the ripples, further 
observed by Ayrton (1910), eroding the ripple troughs and building up the crests. Such 
bedforms, termed ‘vortex ripples’ by Bagnold (1946), exert a much larger drag on the 
flow than friction on sand grains. Vortex ripples occasionally have been ealled ‘orbital 
ripples’ because their wavelength is related to the near-bed orbital diameter of the wave 
motion, or ‘megaripples’ when they exceed some large wavelength, although they should 
not be eonfused with nearshore short-erested megaripples generated by different proeesses 
(Gallagher eta/., 1998). 

Based on dimensional analysis and numerical morphodynamic modeling Andersen 
(1999) and Andersen & Fredspe (1999) found that the wave flow over ripples is essentially 
governed by two nondimensional parameters, X/d and fl/X, where X and T] are the rip¬ 
ple wavelength and height, and d is the diameter of the bottom orbital excursion of water 
pareels (figure 5). Sediment motions are governed by two additional parameters, the ratio 
Ws /Mmax of the Settling velocity of sand grains Ws and the maximum near-bed orbital ve- 
loeity Mmax, and the maximum ratio of frietion and buoyant forees aeting on a sand grain, 
represented by the Shields number \|/max (1-4). When \[tmax is larger than a eritieal value \|/c 
the flow is able to move sand grains. If the bed is initially flat, ‘rolling-grain ripples’ will 
form as a result of an instability of the flat bed. In this ease eaeh grain creates a region 
of weaker flow (‘shadow’) in its lee, and grains tend to group and form ripples with larger 
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Figm e 5. Schematic of ripples 

Definition of ripple and wave forcing parameters. The cartoon conesponds to the passing 
of a wave hough, propagating to the right. Tiubulent eddies are generated in the lee of the 
ripple crests. The wave flow (represented by the fiee-stream velocity //max) and tlie vortex 
retmn flow (dashed anow) converge at the ripple crest, building up the ripple profile. The 
velocity of a water particle (open cucle) at the top of the bomidary layer is indicated by 
//max, and its hajectory is a flat ellipsis. The velocity of a suspended sand gr ain (filled cucle) 
is a combination of the siuToimding flow velocity, comparable to //max, and its settling 
velocity 


shadows (Blondeaux, 1990; \^ttori & Blondeaux, 1990; Andersen, 1999, 2001). These 
rolling-grain ripples evenhrally evolve into vortex ripples (Sherer, Melo & Marder, 1999; 
Faraci & Foti, 2001). The height r| of vortex r ipples is generally closely related to A,. Tire 
ripples ar e steepest for 1 < V)/max/M^c < 4, when the vortices created in the lee of the crests 
maintain \\/X values between 0.1 and 0.15. 

Niunerical model simulations of the morphodyuamics of oue-dunensional ripples 
imder simrsoidal waves confirm that the vortex forming m the lee of the crest with a size 
of the order of the orbital diameter d exerts a strong shear on the lee-side slope of the 
r ipples that tends to build up the crests together with the flow on the upstr eam slope of the 
ripple (Andersen, 1999). If two ripples are initially closer than a minimal stable wavelength 
the vortex in the lee of the upstream ripple will erode the downstream tipple and one 
ripple may disappear creating a default in the regular spacing of the bed that will allow A, 
to grow. Conversely, ripples that initially are much farther apart than Xm will promote the 
generation of ripples in the troirghs, dividing A,by a factor two. Values of A.„ were formd 
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to be related to the orbital diameter with Xm = OAd when Ws/u^ax < 0.07 and = 0.63d 
when Wj/Mmax > 0.07. The transition for w^/Mmax = 0.07corresponds to a shift in the 
principal mode of sediment transport from suspended load to bed load, and the evolution 
from two-dimensional to three-dimensional ripple patterns (Nielsen, 1979). The movement 
of defects in the three-dimensional ripple pattern should tend to reduce the average crest 
distance to Xm (Andersen and Fredspe, 1999). Other classifications and characterizations of 
wave-generated ripples have been proposed based on empirical evidence (e.g. Mogridge, 
Davies & Willis, 1994; Wiberg & Harris, 1994) but they generally failed to reconcile all 
laboratory and field observations. 

In the field ripples may be affected by the presence of wave groups (e.g. Madsen 
etal., 1990), the mixture of grain sizes (e.g. Wallbridge etal., 1999), and the directional 
distribution of the waves (Willis et al, 1993), although experiments support simple param- 
eterizations using ‘equivalent parameters’, e.g. the median grain size D 50 , and the velocity 
Ms, orbital diameter d^ and Shields number \|/s based on the significant wave height (e.g. 
Traykovski et a/., 1999). 

The presence and characteristics of wave-formed ripples on the North Carolina 
continental shelf is verified in chapter IV with simultaneous wave and bottom morphol¬ 
ogy measurements acquired in 1999 during the SHOaling Waves Experiment (SHOWEX). 
These observations were motivated by earlier analysis of wave data (Herbers, Hendrick¬ 
son & O’Reilly, 2000), and numerical modeling studies (chapters II) of strong damping of 
swell propagating across the shelf. The new data presented in chapter IV includes sediment 
samples, repeat sidescan sonar images of the bottom, and three-month-long observations 
of wave frequency-directional spectra. The analyzed ripples direction and wavelength, in 
relation with preceding wave forcing conditions and sediment properties, are reconciled 
with the parameterization of Andersen & Eredspe (1999) and previous observations by 
Traykovski et al.{l999). 
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3. Parameterization 

The representation of sand bedforms in wave models (e.g. Graber & Madsen, 1988; 
Tolman, 1994) usually involve a ‘ripple roughness predictor’ which, based on the wave 
conditions and sediment nature, determines the flow regime, the type of bottom features 
(e.g. Clifton, 1976; Wiberg & Harris, 1995) and their equivalent sand grain roughness 
ki^ (e.g. Grant & Madsen, 1982; Madsen eta/. 1990; Li & Amos 1998). This roughness 
predictor is combined with a hydrodynamic model of the bottom boundary layer flow that 
predicts the corresponding wave energy dissipation. Most hydrodynamic models param¬ 
eterize turbulence with a vertical profile of the eddy viscosity (Kajiura, 1968; Grant & 
Madsen, 1979; Weber, 1991a, 1991/?; see Wiberg, 1995, for a review). The use of a single 
roughness length for spectral waves was validated in laboratory experiments by Mathisen 
and Madsen (1999). The general parameterization of (non-linear) spectral dissipation term 
(the bottom friction ‘source’ term) in terms of the wave spectrum was considered by Weber 
(1991) and simplified for practical application using a narrow spectrum approximation that 
was also used by Madsen (1994). 

Grant & Madsen (1982) provided the first comprehensive parameterization of the 
interaction of waves with a mobile sandy bed (i.e. relic roughness, active ripples and sheet 
flow). The present dissertation considers a later parameterization proposed by Tolman 
(1994) that combines a ripple roughness predictor proposed by Madsen et al. (1990), with 
Grant and Madsen’s (1979) hydrodynamic model, extended to spectral waves by Madsen, 
Poon & Graber (1988) based on a narrow spectrum approximation. For the sheet flow 
regime Tolman used Wilson’s (1989) extrapolation of river flows to oscillatory boundary 
layers. The source term S is expressed as a quasi-linear function of the energy density E 
with a directionally isotropic (negative) growth rate X: 


5(/,0) 

Hf) 


X(/)x£(/,0) 


(I.5a) 

(I.5b) 


where g is the gravity acceleration, h is the local water depth, k is the wavenumber mag¬ 
nitude, and fe is the local dissipation factor representing the ripples or sheet flow effect 
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depending on sediment and wave characteristics. The sediment parameters are a represen¬ 
tative grain size D, specific density 5 = p^/p, where p^ and p are the densities of sediments 
and water respectively, and the critical Shields number for sediment motion V[tc. The wave 
parameters are a representative orbital velocity Mb,ms and a horizontal displacement Mb,ms 
at the top of the bottom boundary layer, evaluated using (11) and (25) in Madsen et aUs 
(1988) model: 


2 _ 

^b,ms 

./k smh [kh) 

(I.6a) 

‘.*b,ms 

/—^E(k)rfk. 

./k smh [kh) 

(I.6b) 


For a linear profile of eddy viscosity. Grant and Madsen (1979) determined the skin friction 
factor f!^, giving the Shields number \\frms = ms/ [<? (■^ ~ 1) ^]’ total friction factor 
fw (ratio of bulk stress and as implicit functions of the grain size D and equivalent 

grain roughness of the bedforms k]\f, respectively: 

2 D or kj^ 


I 


f'worf^ = 


fiv fw 30k Mb,ms 


(I.7a) 

(1.7b) 


ker^ -bkei^ (2^/z^l^ 

where zo/^ is a nondimensional roughness length, K is Von Karman’s constant (k = 0.4 for 
clear water), and ker and kei are the zeroth order Kelvin functions. 

The dissipation factor fg is assumed equal to the total (skin friction and form drag) 
friction factor /w (e.g. Nielsen, 1992), determined by solving iteratively (I.7a,b) The equiv¬ 
alent grain roughness k^ of the bedforms is parameterized as a function of Vc, Mb,ms> 
and Mb,ms- For Vj/rms/tl/c <1.2 (i.e. in the ‘relic roughness’ regime), kjv is taken to be 0.01 
m and fg is limited to a maximum value of 0.30 (Jonsson, 1980). Beyond 1.2v[tc, in the 
‘active ripple/sheet flow’ regime, ki\/ is the sum of a ripple roughness k^and a sheet flow 
roughness ks. Madsen et a/.(1990) gave empirical values of k^for random waves in labora¬ 
tory experiments and Wilson (1989) extrapolated to waves values of k^ measured for river 
flows: 


kr 


^b,ms 


X 1.5 
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Figiue 6. Examples of dissipation factors fg as a fimction of the Shields number V|/njis 
The solid line is Tohnan’s (1994) parameterization for a representative grain size D = 
0.15 mm. and wave period T = 14 s. Tire dashed line shows conespondurg vahtes of fe 
irsing the JONS WAP parameterization. 




(I.8b) 


The first evahration of Tohuan’s (1994) parameterization against field data is per¬ 
formed in § n.E, followed by a systematic evahration rtsing all swell-dominated data sets 
fiom DUCK94 and SHOWEX in chapter V. Slight modifications to Tohnan’s (1994) for- 
mirlation are proposed in chapter V, with empuically calibrated coefircieuts that improve 
the accruacy of swell hindcasts on the North Carolina continental shelf 

Also considered in § II.E and V is the simpler empuical ‘JONSWAP’ parameteri¬ 
zation of bottom dissipation irsed in many operational wave prediction models. It assiunes 
that fe is inversely proportional to //b,nm so that the atteniration coefficient F = g/e?/b,nns/2 
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is constant and the source term is given by: 


5(/,e) = -r 


271/ 


gsm]\{kh) 


Eif.Q) 


(1.9) 


An average value F = 0.038 m^s”^ was inferred from the JONSWAP North Sea experi¬ 
ment (Hasselmann et al, 1973), despite a wide scatter of the attenuation coefficient F (from 
0.0019 to 0.160 m^s^^). This parameterization has encountered some success (Bouws and 
Komen, 1983), and replaced quadratic drag formulations proposed previously (Hasselmann 
and Collins, 1968). The reason for this success over sandy bottoms, in spite of very few 
physical arguments, probably comes from the fact that this parameterization gives dissi¬ 
pation factors fe that decrease as a function of the Shields number (figure 6), following 
the movable-bed model for relic roughness and sheet flow conditions, although it clearly 
misses the amplification of fe in active ripple generation conditions. 


F. NUMERICAL WAVE MODELS 

Most numerical models for the evolution of surface gravity waves across ocean 
basins, marginal seas, and continental shelves that account for non-conservative processes 
(e.g. wave generation by winds, wave breaking, and bottom friction) are based on imple¬ 
mentations of the spectral energy balance (1.2), or a similar wave action balance, with finite 
difference schemes. These models are efficient in deep water applications where large 
spatial and temporal scales of wave evolution allow for relatively coarse grids (e.g. The 
SWAMP group, 1984; Komen et al, 1994). 

In shallow water accurate representation of refraction may require grid resolution 
of the order of 100 m. If the region of interest is small (less than 100 km^), a high reso¬ 
lution Eulerian model is feasible and gives good results (Booij, Ris & Holthuijsen, 1999), 
but the computational cost is presently too large for larger shelf areas, even in a steady state 
formulation. Additionally, finite difference approximations in these models cause numeri¬ 
cal diffusion, artificially spreading wave energy in time, x, and k space, in a way unrelated 
to the physical evolution of a wave spectrum over bottom topography. High-order finite 
difference schemes and piecewise ray methods, using local ray trajectories to estimate the 
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advection terms of (1.2), have been developed to mitigate this effect (Sobey, 1986; Young, 
1988; Lavrenov & Onvlee, 1995; Benoit, Marcos & Becq, 1996). 

Lagrangian wave prediction models based on (1.3) usually assume a source term S 
equal to zero. This approach is suitable for narrow shelf regions where propagation dis¬ 
tances are too short for significant wave generation or decay (O’Reilly & Guza, 1991). 
Lagrangian models avoid the numerical diffusion of finite difference schemes, but the ray 
trajectories are highly sensitive to topography details. The scattering of rays over rough bot¬ 
tom topography causes physical diffusion of wave energy that may broaden wave spectra 
in shallow water. The accurate representation of these fine scale bathymetry effects in a ray 
model requires averaging over a large number of rays, whether the rays be computed from 
initially parallel directions (forward refraction, e.g. Bouws and Battjes, 1982) or from fixed 
points (back-refraction, e.g. O’Reilly & Guza, 1993). Back-refraction models are not based 
on finite area elements, unlike forward refraction and finite-difference schemes in Eule- 
rian models, and thus have different conservation properties. For example finite-difference 
schemes are generally constrained to conserve energy fluxes through the model domain, 
but the energy fluxes obtained through spatial interpolation in a back-refraction model bal¬ 
ance exactly only in the limit of high spatial and wavenumber resolution. Nevertheless 
if a detailed bathymetry is available, a back-refraction ray model with high wavenumber 
resolution gives a potentially more accurate representation of wave propagation than finite 
difference schemes. 

Cavaleri and Malanotte-Rizzoli (1981) included wind input and dissipation source 
terms in a ray model. Their model parameterizes the source terms for each individual wave 
component, and solves the energy balance equation independently for each ray, without any 
coupling. Lavrenov (Personal communication, 2000) used a similar approach but restricted 
the source terms to those that only depend on the energy at the same wavenumber, such 
as the damping of waves by sea ice, so that rays are indeed independent. In chapter II 
we present a new model that includes coupling between rays through a source term that is 
parameterized in terms of the full energy spectrum. 
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II 


THE CREST WAVE MODEL 


tin this chapter, published in a slightly different form in the Journal of Physical 
Oceanography (Ardhuin, Berbers & O’Reilly, 2001), a new numerical model is presented 
for solving the spectral energy balance (1.3). The source term S (k,t) is evaluated at each 
point of a coarse Eulerian grid, and subsequently interpolated from this grid onto ray tra¬ 
jectories. The energy balance (1.3) is integrated along a full spectrum of rays traced back¬ 
ward from each grid point to the model boundary. Spectral components are advected from 
the model boundary along the precomputed rays while being modified by the interpolated 
source terms, until they reach a grid point where all components are combined into a full 
spectrum E (k, t) from which S (k, t) can be evaluated (figure 7). The advection and source 
term computations are performed simultaneously for the entire model domain. This hybrid 
Eulerian-Lagrangian model essentially couples a Eagrangian energy advection scheme with 
an Eulerian source term computation scheme. The formulation of the source term compu¬ 
tations is not constrained in any way by the advection scheme and thus can be adapted from 
existing third-generation models. 

A. NUMERICAL SCHEMES 

The model consists of two parts. Eirst wave rays are traced backwards from fixed 
Eulerian grid points, with positions x,-, to the model boundary. Second, these trajectories 
are used to integrate (1.3) in time, using an ensemble average over a large number of rays. 
Along each ray, arriving at x,- with a wavenumber vector k, we define a Eagrangian energy 
density (t,x) as the energy density ‘upstream’ of x,- at time t, where x is the energy ad¬ 
vection time from the local ray position to the grid point x,- (figure 7). The spectral densities 
E^ are averaged over ensembles of rays within finite bands k^ of the arrival wavenumbers 
k at x/. The full Eulerian energy density spectrum E^ (x,-,k,t) at x, is evaluated by com¬ 
bining the average Eagrangian density predictions E^ (t, 0) at x,- for all bands k^. A source 
term S (x,-,kj,t) is determined at each grid point from the full Eulerian spectrum E^and 
other local parameters (e.g. wind stress and bottom roughness). S is then interpolated in 
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Figure 7. Numerical scheme of CREST 

The Lagiaugian energy balance is integrated from t to t+Af along a single ray (solid ciu've) 
using a spatially inteipolated somce teim. Filled cucles symbolize the magnitude of the 
energy density, and dashed arrows mdicate the interpolation of the somce term fr om the 
Eulerian grid (sqirares) onto the ray at increments 6x. See § II.A for frrrther details. 


X and k space to yield an approximate somce term 5(t,T) at the local ray positions and 
waverrmnbers which in hmi modifres along the rays (frgme 7). Rays and grid are 

thus coitpled at x = 0 only. 
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The entire set of interpolation coefficients, representing the influence of the topog¬ 
raphy on waves is precomputed once, and stored in files. Using these files and a time series 
of wave spectra at the model open boundaries the energy balance equation is integrated in 
time. 

Although the Lagrangian energy balance (1.3) holds only for energy density in 
wavenumber (k) space, the propagation of waves is formulated more conveniently using 
wave frequency / and direction 0 as variables. In the following / and 0 are used through¬ 
out in ray calculations, grid discretization and result displays, but the energy density in k 
space is used in the energy balance calculations. 

1. Model domain and boundary conditions 

The model domain covers a region of known bottom topography. From an arbitrary 
set of Ngp grid points (hereinafter called ‘model grid’) with locations (x,);^j a triangular 
mesh is generated using Delauney’s tessellation technique. The outermost points of the 
mesh form the model boundary, which is therefore a polygon. Additional interior polygons 
can be added to the boundary in order to represent islands in the model domain (figure 
8). Ray trajectories are traced backward in time from the grid points X; until they cross a 
boundary. For each x,, rays are computed for a large number of frequencies fj and arrival 
directions 0/. Depending on the geographical region covered by the model domain, rays can 
be trapped in shallow water and end at the coast, or reach a deep water region where they 
become straight, or cross the model boundary in a region of intermediate depth. In all cases 
a ray is terminated when it crosses a triangle side connecting two boundary grid points, and 
the Lagrangian energy density carried by the ray into the model domain is approximated 
by a linear interpolation of the spectral densities at these two grid points (figure 8). 

The boundary condition for the model is therefore fully prescribed by the spec¬ 
tral densities at the grid points along the boundary, for directions toward the inside of the 
boundary. On the open part of the boundary, spectra may be estimated from deep water 
wave measurements or obtained by nesting the model within a larger scale wave model. 
On the closed coastal part of the boundary, the energy entering the domain may be set 
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Figiue 8. Treatment of the boimdaiy condition 
In this schematic squares represent grid points. The boimdaiies (dashed lines) separate and 
couple adjacent model subdomains. Examples are shown of rays transporting energy into 
the model domam from foiu different types of boimdaiies (a: shelf break, b: island, c; coast, 
d: internal boimdaiy between model subdomains) to a given grid point (large square). In 
cases a and d energy is advected through the boimdary, whereas in cases b and c energy is 
reflected fr om the boimdaiy. hi all cases, the energy is inteipolated (dotted arrows) at the 
boimdaiy from the adjacent two boimdaiy grid points. 
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equal to zero (i.e. wave energy impinging on the coast is dissipated in the surf zone) or, 
in the case of a steep coastline, determined by partially reflecting the shoreward energy 
flux. In order to reduce the scattering of rays over large propagation distances, the model 
domain can be subdivided into subdomains that are coupled through their common bound¬ 
aries. This technique reduces memory requirements by shortening the rays, at the expense 
of some local numerical diffusion, as the energy that is transmitted through the boundary 
is interpolated from boundary grid points (e.g. ray d in figure 8). 

2. Precomputations 

a. Rays 

In applications presented here the model domain is small enough to neglect 
the curvature of the earth, and use local Cartesian (.r,y) coordinates. The geometry of wave 
rays is determined by Fermat’s geometrical optics principle that the integral of the phase 
speed C along a curve is minimum when this curve is a ray, which yields Snel’s law ^ when 
bottom contours are parallel. The ray equations are: 


dx 

ds 

= cos(0) 


(II. la) 

dy 

ds 

= sin(0) 



(II. lb) 

dd 

ds 

1 dC 

C dh 

dh . . , dh 
— ■ sm(0) —— - cos (0) 
dx dy 


(II. Ic) 


with s a curvilinear coordinate along the ray, h the water depth, and 0 the angle between the 
X axis and the tangent to the ray. Wave energy is transported along the ray with the group 
velocity Cg and the frequency / is conserved. In the linear approximation we have: 

{2nf)^ = gA:tanh (kh) (II-2) 


'The dutch mathematician Willebrord Snel discovered the law of refraction in 1621, but it was only 
published in 1703 in Dioptrica, by Christiaan Huygens, in which Snell is given the Latin name Snellius, 
which is often misspelt in English as ‘Snell’. The French philosopher Rene Descartes gave Snel’s law in La 
dioptrique, an appendix to his famous Discours de la methode pour Men conduire sa raison et chercher la 
verite dans les sciences, published in Leiden in 1637, but it was apparently taken from Snel’s work although 
he repeated Snel’s experiments in 1626 or 1627. Snel’s law is wrongly attributed to Descartes in the French 
scientific literature. 

Source: the MacTutor history of mathematics archive, University of St Andrews, Scotland, http://www- 
groups.dcs.st-andrews.ac.uk/ history 
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c 


Q = 

where k= |k| is the wavenumber magnitude. 

Along the ray the local depth and bottom slopes are evaluated from a bi¬ 
quadric fit to the bathymetry grid (Dobson, 1967). The wavenumber magnitude k is com¬ 
puted from / using (II.3) and used to determine C, Cg and With these parameters (II.I) 
are integrated using an error-controlled Cash-Karp Runge-Kutta scheme (Press et al, 1992) 
with a variable step size. 

Along each ray the position and direction (x™, 0“) are computed at small 
distance intervals 55 = ^^^Cgdx that correspond to a fixed advection time step 5x. A 
5x was chosen for each frequency such that 5^ = 200 m in deep water. The result of the 
ray computation is a series of positions and directions (x™, 0“) for each of the rays with m 
ranging from 0 at the initial grid point to M at the domain boundary, with typical values 
M ~ 1000 in the implementation presented in § II.D. 

M, and 0^ give the time lag, position and direction at the end of the ray, 
needed to specify the boundary condition. Although waves can travel along the same ray 
in both directions, the rays are used here only to advect energy from the boundaries to the 
grid points. 

b. Interpolation of boundary conditions and source term 

At each position x™ along a ray, a local source term estimate S (x™) is given 
by the linear interpolation in space of source term predictions at the three grid points x,- of 
the local triangle. Since the source term is computed only at discrete directions 0/, another 
linear interpolation, with weights wf, is performed over the two directions 0/ that enclose 

the local direction O'” of the ray. The same procedure is used for deriving an estimate Eg 

of the energy density (x^) at the boundary: 

5(x"') = Y.a’fwfSixM (II.5) 

ij 

Eb = £p,wf£^(x,-,0;) (II.6) 
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where the spatial weighting coefficients a™ and (3/ are nonzero only for the three grid points 
X,- on the vertices of the local triangle, and the two points of the boundary segment crossed 
by the ray, respectively. 

In order to resolve the refraction of a single wave component and interpolate 
accurately the source term onto the ray, a small time step 5x is required that is of the order 
of 10 to 100 s for typical swell group velocities (Cg = 0(10 m s^^)) and scales (0(1 — 
10 km)) of bottom features. This time step is too small for an efficient time integration of 
the energy balance equation (1.3). This integration is performed here with a fixed larger 
time step At (10 minutes in the calculations presented here), that resolves the typically 
slower evolution of the wave energy and source terms in space, and temporal changes of 
the offshore boundary conditions. The source term S (II.5) is averaged over an advection 
time interval At”, that covers values of x from (n — 1) At to nAt: 

S" = Y.A^iS{xM (II.7) 

id 

5t 

with A"; = ^ 

i'=i,l'=l,m 

where the summation over m includes all ray segments that fall within the time step At”. 
In the applications presented here the timestep index n ranges from 0 at the grid point 
to 10-50, depending on the location of the grid point, the frequency of the waves, and the 
complexity of the topography. Higher frequency waves and rough topography require more 
timesteps than low frequency waves and smooth topography, because the group velocity 
decreases with increasing frequency and bending of rays over rough topography lengthens 
the propagation path. 

c. Finite bandwidth approximation 

So far we have considered the evolution of the spectral energy density {t, x) 
along a single ray. Since individual ray trajectories are highly sensitive to the underlying 
bathymetry, the energy balance equation is ensemble-averaged over a ‘bundle’ of rays orig¬ 
inating from Xi with frequencies and directions covering a small but finite bandwidth. The 
rays that form a bundle can be scattered and follow different paths away from x,-, there¬ 
fore the ray ensemble has a physical interpretation only at the grid points x/ as a finite 
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bandwidth average. The ensemble averaged energy density and associated source term 
interpolation coefficients A"; are given by weighed averages of single-ray values: 

= '^brE^it.x) (IL9) 

r 

A-i = I Mr (11.10) 

r 

where the summation is over all the rays in the bundle, and br is the fraction of the finite 
bandwidth attributed to the individual ray r. 

Different rays from the same bundle may reach the boundary during differ¬ 
ent timesteps, so that the ensemble average ‘boundary energy’ Eg must be defined for each 
time step n\ 

Ei = ( 11 . 11 ) 

ij 

B'l, = E (11.12) 

r,i'=i,l'=l 

where the summation is restricted to those rays that reach the boundary during time step n. 

Averaging over finite frequency-direction bands not only accounts for the 
scattering of rays by refraction over bottom irregularities, but also has the advantage of 
avoiding the ‘garden sprinkler effect’ of Eulerian models formulated for a discrete spectrum 
(e.g. SWAMP group, 1984). A large number of ray computations (of the order of 1000 for 
applications presented here) may be needed to obtain a stable ensemble average but these 
time-consuming computations can be performed in parallel for different bundles and grid 
points. 

The results of the precomputation are the ensemble averaged interpolation 
coefficients and B^. Theses coefficients are written to files that are used in the time- 
integration scheme described below. 

3. Integration in time 

The energy balance equation (1.3), averaged over ray ensembles, is a unidimen¬ 
sional time evolution equation that can be integrated using standard finite difference schemes 
However, a more accurate formulation is possible for linear or quasi-linear source terms. 
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The total source term S is split into a (quasi-)linear part XE and a residual term R that 
includes constant and non-linear (in E) contributions: 

— = IE + R (11.13) 

at 


For X and R constants, (11.13) has an exact solution for the evolution of E over one time 
step: 


E (t At) 


E (t)exp(?iAt) +R 


exp (XAt) — 1 
X 


(11.14) 


For X and R varying slowly in time and space an approximate solution is obtained by 
replacing X and R in (11.14) with average values. The interpolation of the total source term 
XE-\-R h more accurate with this formulation, provided that the gradients of X, in k-space, 
x-space and time are smaller than those of E (see § II.C). For fully non-linear source terms 
(i.e. X = 0) (11.14) reduces to a first order Euler scheme. 

X and R are assumed to be known functions Q. and 0 of the local wave spectrum 
that can be adopted from parameterizations in existing Eulerian models. X and R are inter¬ 
polated from the Eulerian grid onto the rays using the precomputed coefficients The 
complete integration scheme is given by: 


Source term evaluation (on the grid) : 

X{t) = a{E^ (t)) 

(11.15a) 

R(t)=0(E^ (t)) 

(II. 15b) 

Interpolations (grid to rays coupling) : 


£SW = Es';£«(x,-,e,,() 

i,l 

(II. 15c) 

r{t)= Y.A^M^iME^{xiA,t) 

(II.15d) 
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Prognostic equation (along the rays) : 


(Il.lSe) 



+R"it) 


exp (?) Atj — 1 


(IL15f) 


Rays to grid coupling (at x = 0) : 

(? + A?,0) 


E^{t + M) = { 


or 

Eb (t + A?) 


(II.15g) 


where the frequency variable fj is omitted. Variables x/, and 0/ are written explicitly only 
in the interpolations (IL15c-e). In (IL15d) the weighting of X by the corresponding energy 
density E^ allows the conservation of the source term XE in the interpolation. The prog¬ 
nostic equation (II.15f) applies the interpolated boundary condition and source term to the 
Lagrangian energy balance to determine E^ at the next time step. The Eulerian spectrum 
E^ is advanced to time ? -I- A? with (II.15g), closing the set of equations. For grid points x, 
located on model domain boundaries, the spectral densities E^, for waves traveling into the 
model domain, are prescribed by the boundary condition Eg. On the deep water boundary 
Eg is set equal to the observed deep water spectrum. At other external boundaries Eg is set 
equal to zero. At internal boundaries E^, for waves traveling into one domain, is prescribed 
by E^, for waves traveling out of the other domain. For all other components and interior 
grid points E^ follows from E^. Each equation can be evaluated in parallel for all the ray 
ensembles and all grid points, and different frequency bands are only coupled by the source 
term. 


The accuracy of this scheme depends on the relative size of the Eulerian {Te) and 
Eagrangian {Tg) time scales of wave evolution. For Tg Tg (e.g. a sudden and uniform 
change in forcing conditions over the entire model grid) the dominant source of error is the 
low order time integration scheme. If Tg 3> Tg (e.g. strong energy dissipation at a fixed 
location, with quasi-stationary boundary conditions and source term) the largest errors may 
result from spatial interpolation of the source term. Earge errors occur if either Tg or Tg are 
comparable to, or smaller than. A?. For all cases presented in § ILF, A? is small compared 
to both Tg and Tg. An alternative predictor-corrector scheme was tested, giving results that 
are indistinguishable from those of the scheme used here. 
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B. SOURCE TERM 

In the present chapter the source term is restricted to the energy dissipation caused 
by bottom friction. The detailed parameterization, that represents the rate of change of the 
energy of the spectral components, is taken from Tolman (1994), as presented in chapter I 
and discussed again in chapter V. The numerical scheme of the present model assumes that 
this source term varies smoothly on the scale of the Eulerian grid. This assumption may be 
violated in the case of bottom friction, and the spatial interpolation of the source term onto 
ray trajectories may cause large errors in the transition region from the ‘relic roughness’ to 
the ‘active ripple’ regime where fg increases by one order of magnitude (figure 6). Tolman 
(1995) proposed a subgrid model of the source term that accounts for subgrid variations 
of the Shields number \\f resulting from variations over each grid cell of the water depth h, 
grain size D, critical Shields number for sediment motion \|/c, significant wave height H^, 
and peak wave period Tp. For simplicity these five random variables were assumed to be 
Gaussian and independent. Because no quantitative information on the spatial variability 
of sediment characteristics was available, a simpler subgrid model was implemented here. 
If uniform sediment properties are assumed then both D and \|/c are uniform within each 
grid cell, leaving only three random variables h, and Tp. In model simulations of swell 
evolution on the North Carolina shelf, using observed incident wave conditions, most of 
the subgrid variability of the source term resulted from the subgrid variations of water 
depths h rather than the wave parameters and Tp, and was correlated with h. This 
predominance of the depth variability was also noted by Tolman (1995) and used in his 
computations. 

In the present subgrid model spatial variations in the water depth h are represented 
by forming a histogram of depths for a grid cell that consists of the triangles surrounding 
the grid point x/, using ten depth bins that span the mean ± 2 standard deviations of h. A 
corresponding linear theory shoaling correction of the wave spectrum is added to account 
for correlations of with h. A subgrid-averaged value of X is obtained by averaging 
estimates of X (1.5b) for each depth bin, based on the corresponding shoaling-corrected 
wave spectrum, weighted by the depth histogram values. 
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There is no strong justification for neglecting the subgrid variability of the sediment 
properties and simplifying sediment description by using the median grain size and a 
constant critical Shields number \|/c = 0.05. This problem is further discussed in chapter IV 
where we show evidence of significant changes in sediment properties over short distances 
(O (100 m)) on the North Carolina continental shelf. 

C. MODEL TESTS 

The error in the ray computation is controlled by the variable time step, but other 
errors are introduced by the discretization in frequency and direction and the ray ensemble 
average. The accuracy of the propagation scheme was tested by applying the model, with 
the source term set equal to zero, to an idealized shelf with parallel depth contours, for a 
uniform and stationary offshore boundary condition (figure 9a). The mean wave directions 
and directional distributions of energy predicted by the model agree closely with analytical 
(Snel’s law) results (figure 9 a,c), demonstrating that ray integration and discretization 
errors are small. 

The model formulation assumes that the source term varies on scales comparable 
to or larger than the spacing of the grid. This condition is required for a valid interpolation. 
The accuracy of the interpolation scheme was tested by computing the source term directly 
at 20 additional grid points located along a ray segment (for 0.07 Hz waves) that covers a 
full source term integration time step At = 600 s. Source term estimates interpolated onto 
this ray segment with (II.15d,e) in a hindcast of wave evolution across the North Carolina 
shelf (§ II.E), are compared with direct estimates at the additional grid points in figure 
10. Results (averaged over a time step At) show that the linear spatial interpolations give 
a good approximation of subgrid variations in the source term. That is, the source term 
gradients are rather well resolved by the grid. The interpolation is most accurate when 
the ripple regime is the same at all the neighboring grid points. Overall, a quasi-linear 
implementation of the movable bed source term (i.e. X = S/E, R = 0 in (11.13)) (figure 9a), 
yields smaller errors than a non-linear implementation of the same source term (i.e. X = 0, 
R = Sin (11.13)) (figure 9b). 


36 






North 30 60 East 120 150 South 

directkxi (from) 


Figuie 9. Model tests with an alongshore imifoiiu shelf 
The offshore wave energy is distributed uniformly over a naiiow (0.0655-0.0685 Hz) fre¬ 
quency band, (a) Cross-shelf depth profile, (b) Predicted (+) and analytical (solid line) 
mean wave direction versus cross-shelf distance. Results are shown for naiiow offshore di¬ 
rectional distributions (standard deviation of 10 degrees) with mean wave directions vary¬ 
ing between 30° and 150°. (c) An example directional distribution predicted by the model 
(+) in 20 m depth for a given offshore bimodal distribution (dotted line) is compared to the 
analytical solution (solid line). 
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S/E : Direct computation (1/hour) 


Figure 10. Comparison of semi-implicit and explicit schemes 
Interpolated versus directly computed values of S/E (inverse of the bottom dissipation e- 
folding time), using (a) a quasi-linear source term implementation that interpolates X (L5b), 
and (b) a non-linear implementation of the same source term that interpolates S (L5a). The 
source term estimates are averages over a 10 minute time step. Symbols represent different 
boundary layer regimes within the grid cell: relic ripples (-t), active ripples (A), and a 
transition from relic to active ripples (squares). 
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D. MODEL IMPLEMENTATION AND EIELD DATA 

The hybrid Eulerian-Lagrangian model ( § II.A) with a movable bed bottom dissi¬ 
pation source term ( § II.B) was implemented for the North Carolina shelf region between 
35° N and 37° N (figure 11). During the DUCK94 experiment a 100 km cross shelf tran¬ 
sect of nine bottom mounted pressure sensors was deployed extending from 12 m depth 
(site A) to 87 m depth (site I) (figure 11; Herbers et al, 2000). The instrument deployed at 
site H in 49 m depth was located within 2 km of 3-m discus buoy 44014 operated by the 
National Data Buoy Center (NDBC). Between site A and the shore a pressure sensor array 
was operated in 8 m depth by the Army Corps of Engineers Eield Research Eacility, Duck, 
North Carolina. Other instruments on the inner shelf included current meters, thermistors 
and conductivity sensors in depths ranging from 4 m to 26 m. Data from these instruments 
show that outside the surf zone the depth-averaged currents are mostly wind-driven with 
speeds usually in the range 10 to 20 cm s^^ and occasional stronger currents (> 50 cm s^^) 
in storm conditions (Eentz etal., 1999). These currents are generally much weaker than 
both the wave speeds and the near-bed orbital velocities of energetic swell in shallow water. 
The effects of currents on the propagation of swell and on the wave bottom boundary layer 
are neglected here. 

Bathymetry data was derived from the National Ocean Service digital database and 
additional bathymetric surveys conducted during the DUCK94 experiment (Herbers et al, 
2000). These data sets were merged and linearly interpolated onto a regular 6” longitude 
by 6” latitude grid, using a standard Delauney tessellation technique (11). This grid was 
then linearly transformed into a Cartesian x (west-east) and y (south-north) coordinate grid 
(resolution 152 and 185 m respectively) that is used in the Eagrangian ray integrations (eqs. 
II. 1). Errors introduced by neglecting the curvature of the earth, are small for the size (128 
by 211 km) of our domain (O’Reilly & Guza, 1993). 

The Eulerian model grid, much coarser than the bathymetry grid, consists of 329 
grid points arranged in triangles that vary in size from 5 km on the inner shelf to 10 km 
on the outer shelf (figure 11). Slightly coarser grids gave similar results, suggesting that 
the resolution chosen here is adequate. The model domain was made as small as possible 
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Figiue 11. Bottom topogiaphy and model giid 
The giid points where the soiuce temi is evaluated are the nodes of the triangular- mesh. 
A linear interpolation is applied in each triangle to approximate the somce temi along 
the rays (figiue 1). Tire entire model domain is subdivided into subdomains separated by 
thicker lines. Grid points denoted with dots labeled A to I are the locations of pressiue 
sensors deployed diuing the DUCK94 field experiment. 
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while covering the shelf region through which most of the wave energy, measured by the 
pressure sensors, has propagated. The model boundaries were chosen to be the 11m iso¬ 
bath (except around the 8 m array, where a grid point is collocated with the array), the 400 
m isobath, and the 35^12’ N and 37"58’ N parallels. Swell energy enters the model domain 
only through the deep water boundary where the spectrum Eg is assumed to be 

spatially uniform. The model domain is subdivided into a main subdomain, around sites 
A to I and additional subdomains (figure 11). This allows for the representation of waves 
coming from high incidence angles, and reduces the memory required to store all the in¬ 
terpolation coefficients, including those for trapped rays, to one gigabyte. Trapped rays 
are not necessary for the applications presented in § II.E, since energy enters the model 
domain only through the deep water boundary, but were implemented for applications with 
other source terms (chapters III and V) and reflective boundaries. The use of subdomains, 
described in § II.A, introduces some numerical diffusion for waves propagating across the 
internal boundaries, but these waves, with high incidence angles, generally carry a small 
fraction of the total energy in the present applications. 

For each grid point x,, rays are initially traced for 162 frequencies, at arrival di¬ 
rection intervals of 0.25°. For each frequency the arrival directions are subsequently bi¬ 
sected (O’Reilly & Guza, 1993), until neighboring rays have directions and positions at the 
boundary within 2° and 5 km of one another, respectively. If the number of rays for a 3° 
sector exceeds 500 the bisection is stopped. The rays are grouped in 19 x 120 bundles, that 
represent finite bandwidths of the spectrum (x,-,/j,0/,t) with 19 frequencies fj spaced 
exponentially with a 5% increment from 0.05 Hz to 0.12 Hz, and 120 directions 0/ spaced 
linearly over a full circle with a 3 degree resolution. The number of rays per bundle varies 
from Nf X 13 (initial number of rays before bisecting) to Nf x 500 (an upper limit set for 
broadly scattered bundles), where Nf is the number of frequencies per frequency band. Nf 
decreases from 9 for 0.05 Hz to 3 for 0.12 Hz. 

Wave frequency spectra E integrated over directions, are estimated from 

the measured bottom pressure records at sites A-I, using a linear theory depth correction. 
The offshore frequency-direction spectrum Eg 0/, t) is estimated by combining the fre- 
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quency spectrum obtained from the pressure sensor at site H, with directional distributions 
estimated from the nearby NDBC buoy cross-spectra using the Maximum Entropy Method 
(Lygre & Krogstad, 1986). This spectrum is back-refracted from site H to deep water, 
assuming parallel bottom contours, and neglecting the time lag between site H and the 
offshore model boundary. Although the offshore conditions generally varied slowly on 
time scales of several hours, this spectrum is determined at 10 minute intervals, in order to 
match the model time step At. Frequency-directional wave spectra on the inner shelf were 
estimated from the 8 m depth array near site A (Herbers, Elgar & Guza, 1999). 

Based on core samples collected in 1997 on the inner shelf (Rebecca Beavers, Duke 
University, personal communication, 1999), and earlier geological data covering the entire 
shelf (Milliman etal., 1972; Swift and Sears, 1974) we crudely approximate the bottom 
sediments in the entire area encompassed by the model with a thick uniform layer of fine 
quartz sand {s = 2.65), with a representative grain size D = 0.15 mm and a critical Shields 
number \|/c = 0.05. These approximations are further discussed in chapter IV. 

E. HINDCASTS 

Hindcasts are presented for two time periods in 1994, October 17 through 21 and 
November 16 through 19, that are representative of fall weather patterns causing large 
waves on the North Carolina coast. Wind sea and swell were observed in October, gener¬ 
ated by a storm that moved across the eastern United States into the North Atlantic, whereas 
in November large swells arrived from Hurricane Gordon which remained in the western 
Atlantic, south of Cape Hatteras (see figure 43 in chapter V). In both cases, the model was 
run both with and without the bottom dissipation source term. Runs without the source 
term isolate the effects of refraction and shoaling in the evolution of wave spectra, and the 
difference between runs with and without the source term can be used to assess energy 
dissipation caused by bottom friction. 

1. October storm 

On October 14 and 15 local winds from the north-east were strong enough to con¬ 
tribute significantly to the energy balance on the shelf at the dominant (8-10 s) wave periods 
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(figure 12). Strong wind forcing is evident at the NDBC buoy where the mean wave direc¬ 
tion 0 (defined here as the direction of the first-order moment vector / Jk/A: £ (k) dk^dky) 
follows the local wind direction. As wave generation is not represented in the present 
model, predictions are not expected to be accurate during the spin-up of this storm. 

On October 15 the significant wave height observed at site H reached a maximum 
value //s = 5.3 m (4.3 m in the restricted model frequency range), with a peak period 
Tp = 11 s (figure \2b). After October 16 local winds subsided and Tp shifted to 15 s, 
indicating a transition from wind sea to swell. Between October 15 and 18, //s decreased 
to 2.3 m (time I) followed by an increase to 2.8 m on October 19 (time II), with a narrow, 
swell dominated spectrum (not shown). After October 19 Hs and Tp gradually decreased to 
0.6 m and 10 s respectively. 

Model predictions are presented only for the swell-dominated period October 17 
to 21. Predictions of with bottom dissipation are generally in good agreement with 
observed on the outer (e.g. figure 13a) and inner shelf (e.g. figure \3b). The model 
predicts the expected turning of 0 towards the shore-normal direction, owing to refraction 
by the large scale shelf slope (figure 14). The observed shift in 0, up to 25 degrees between 
the offshore buoy and the nearshore (8 m depth) array, is reproduced by the model (figure 
13c). Observed and predicted 0 in 8 m depth differ by less than 5 degrees. 

Model predictions without bottom dissipation show a small decrease in wave height 
across the shelf that is caused by refraction and shoaling effects (figure \3a,b). The model 
with movable bed friction predicts a strong attenuation of across the shelf (figures \3b, 
14) that is comparable to the observed attenuation. The observed and predicted decay 
across the outer shelf is negligible except for a slight (10%) decrease of on October 19 
when swell energy was maximum (time II in figure 13a). Strong decay of is observed 
and predicted across the inner shelf, up to a factor 2 (equivalent to a 75% energy reduction) 
(figure \3b). The observed and predicted decrease of across the shelf is generally smaller 

when T/s is smaller (e.g. compare times I and II in figures \3a,b). On October 21 when 
was less than 0.5 m, the observed and predicted (with and without bottom dissipation) 
are nearly uniform across the shelf. 
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Figure 12. Wind and wave conditions, October 1994 
Three-hour averages, (a) Wind speed (solid line), wind direction (squares), and mean wave 
direction (x) measured by NDBC buoy 44014 (near site H) (b) Significant wave height and 
peak period estimated from pressure sensor H (peak periods were replaced by the NDBC 
buoy values when smaller than 8 s). Vertical dash-dotted lines labeled T’ (17 October at 
23:30 EST) and ‘IF (19 October at 08:30 EST) indicate times for which more detailed 
analyses are presented in figures 14, 15. 
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Figure 13. Hindcasts, October 1994 

Three-hour averages of observed (solid line) and predicted (-t: with movable-bed source 
term, O: without) significant wave heights at sites F (a) and B (b) after the October storm. 
The dotted lines represent model results at the same sites based on the JONSWAP linear 
damping formulation. The offshore is indicated with a dashed line, (c) The mean 
wave direction 0 measured at the 8 m depth array (solid line) is compared to the model 
prediction with the movable bed source term (-t). The offshore 0 is indicated with a dashed 
line. Vertical dash-dotted lines labeled ‘F and TF indicate times for which more detailed 
analyses are presented in figures 14, 15. 
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Figiue 14. Hiiidcasts ofHs and 0, October 19, 1994 
Model predictions (with movable bed soiuce term) of significant wave height (colors) and 
mean wave direction (arrows) at time II (figiues 12, 13). Dotted lines indicate the 30 m and 
50 rn isobaths. 
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Details of the representation of bottom friction in the model are illustrated in figure 
15. The predicted variation of the dissipation factor fg on the scale of the grid resolution 
confirms the importance of subgrid modeling of the movable bed (Tolman, 1995). Predicted 
boundary layer regimes (relic roughness or active ripples) are sensitive to the offshore wave 
conditions. On October 17 (time I in figures 12, 13), the model predicts relic roughness 
over most of the shelf with dissipation factors fg close to the relic regime minimum (/g = 
0.04) and a sharp transition to active ripples (0.08 < /g < 0.18) in depths shallower than 
25 m (figure 15a). The boundary between active and relic ripples generally follows the 
depth contours. The corresponding local decay rate |?i| proportional to fgUi,sm\r^{kh) 
(I.5b) is enhanced not only by the large increase in fg, but also by the increase of the 
bottom orbital velocity Ub in shallow water. Seaward of site D, predicted with and 
without bottom dissipation are nearly equal to the observed whereas further inshore, 
predicted with and without bottom dissipation diverge sharply and predictions of 
with bottom dissipation reproduce the observed decay of across the inner shelf (figure 
I5b). The strongly enhanced dissipation predicted by the movable bed model on the inner 
shelf is consistent with the observed variations in wave heights. However, the JONSWAP 
parameterization also gives reasonable predictions of in this case. 

On October 19 when the swell energy was maximum (time II in figures 12, 13), 
the movable bed model predicts active ripples on the entire shelf (figure (figure 15c). The 
corresponding values of fg are maximum close to the shelf break (0.1 < fg < 0.12), and 
decrease inshore (fg = 0.04 at site B). A strong decay of wave energy inshore of site G is 
evident in the difference between model predictions of with and without bottom dissi¬ 
pation and these energy losses (on average 0.35 W m ^ over the entire shelf) are consistent 
with the observations (figure 15J). Inshore of site D the model with the source term under¬ 
predicts Hs (overpredicts decay) by about 25 to 50 cm. The JONSWAP parameterization 
on the other hand overpredicts by 40 to 60 cm, as might be expected from figure 6. 
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Figiue 15. Ripples and dissipation factors, October 1994 
(a) Ripple regime based on local mean water depth (hatched for active ripples, blank for 
relic roughness), and dissipation factor fe (contoin interval is 0.02) at time ‘F (figures 
12, 13). (b) Obsei'ved (solid line) and predicted (O: without bottom dissipation, +: with 
movable-bed soiuce temi, dotted: with JONSWAP somce term) at time I, as a fimction 
of cross-shelf distance, (c) Same as (a) for time II (figures 12, 13). (d) Same as (b) for time 

n. 
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2. Hurricane Gordon 

Although the eye of Hurricane Gordon remained south of Cape Hatteras, local 
winds were strong (10-15 ms~^) on November 17 through the morning of November 18 
(figure 16a). During the peak of this event (time III) when Hs — 10 m and Tp ~ 15 s (fig¬ 
ure 16a), the local wind speed was about 13 m/s and the mean wave and wind direction 
differed by about 30 degrees. Estimated values of the wind energy input in the model fre¬ 
quency band (The WAMDI Group, 1988, (2.9)) are below 2 W m^^ on most of the shelf, 
while the predicted bottom dissipation rate is generally between 2 and 10 W m^^ (both 
terms are maximum near the coast). Hence, although bottom dissipation appears to be the 
dominant source term, neglecting the wind input in this case may cause significant errors. 

At time III the model predicts a gradual turning of the mean wave direction from 
115° in deep water to 88° in 8 m depth (figure 17), in good agreement with the mean 
direction (88°) observed at the 8 m array (not shown). Model predictions without bottom 
dissipation yield a decrease in from 8.5 m at site I near the shelf break to 7.4 m at 
site B on the inner shelf. This attenuation, associated with the time evolution of the storm 
and conservative shoaling and refraction processes, accounts for only part of the observed 
decrease of T/s to 5.8 m at site B. Including movable bed dissipation brings the model in 
better agreement with the observations (hgure 18/?). The predicted values of the dissipation 
factor fg are about 10 times smaller than the values for the October event, owing to larger 
Shields numbers. On most of the shelf, fg predictions vary between 0.01 and 0.02 (hgure 
18a), corresponding to sheet how. Active ripples are predicted close to the shelf break in 
depths greater than 40 m. The representative bottom velocity (a linear function of for 
a given spectral shape) is 3 times larger than in the October event. In an absolute sense, the 
dissipation rate |5| 1.5a, proportional to fgu\, is a factor 3 larger than in the October event, 
but the relative decay rate |/i| (I.5b), proportional to fgUb, is a factor 3 smaller. As a result, 
the predicted relative decay of across the shelf, due to bottom dissipation, is weaker 
for the Hurricane Gordon case than the October swell event (a 14% decrease compared 
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Figure 16. Wind and wave conditions, November 1994 
Same format as figure 12). The vertical dash-dotted line labeled TIT (18 November at 
08:30 EST) indicates the time when was maximum offshore. 
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Figiue 17. Hindcasts of i/s and 0, November 18, 1994 
Model predictions (with moveable bed soiuce term) of significant wave height (colors) and 
mean wave direction (anows) at the peak of Himicane Gordon (time ‘IIP in figrue 16). 
Dotted lines indicate the 30 m and 50 m isobaths. 
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Figme 18. Ripples and dissipation factors, November 1994 
Time TIP in figiue 16. Same foimat as figme 15, with the addition of the fe — 0.01 contom, 
and horizontal hatches for sheet flow conditions deteiinined fiom the threshold criterion 
= 0.1where D is m cm (Li and Amos, 1999). 
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to 36% in October, cf. figures 15J and ISb). The JONSWAP parameterization yields 
predictions for Hurricane Gordon, that are close to both observed and movable bed 
predictions (figure 18h). 

F. DISCUSSION 

I. Movable bed model 

The comparisons between observations and model predictions suggest that the ob¬ 
served decay of swell energy across the shelf is primarily the result of refraction and energy 
dissipation in the boundary layer over a movable sandy bed. Predicted wave frequency 
spectra (not shown) are also in good agreement with observed spectra, except at very low 
frequencies (/ < 0.06 Hz) where energy levels are relatively low. The hindcast results 
suggest large spatial and temporal variations of the dissipation factor fe as the seabed tran¬ 
sitions through different roughness regimes (figure 6). Tests with different sand grain sizes 
in the range of probable values for the region (0.15 to 0.2 mm) indicated little sensitivity of 
the results. Although more accurate offshore wave data and detailed sediment distributions 
are needed for comprehensive tests of the bed roughness parameterization, the present re¬ 
sults show a model tendency to overpredict swell damping, in particular in the active ripple 
regime (figures 13, 15J). The parameterization of the ripple roughness hr was tuned to 
reproduce laboratory experiments with irregular but unidirectional waves (Madsen et ah, 
1990). Field conditions, with directionally spread waves, are likely to generate more ir¬ 
regular and less steep ripples, with smaller roughness kr, than laboratory experiments (e.g. 
Nielsen, 1981). Therefore the estimates of kr may be biased high. A reduction of ky by 
30% significantly improved the model accuracy (not shown). 

The JONSWAP parameterization gives a relative decay in wave height 
^/^offshore that is Constant for a given dominant frequency. The observations presented 
here, all for swell with a peak period Tp ~ 15 5, instead show a variable relative decay, in 
response to changes in the wave height (e.g. figure 13h). Equivalent values of F, the JON¬ 
SWAP coefficient, were inferred from the movable-bed model hindcasts. At site F on the 
outer shelf, we find 0.025 m^s"^, 0.11 m^s^^, and 0.050 m^s~^ at times I, II and III respec- 
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lively. For the same times at site B on the inner shelf, F values are 0.11 m^s“^, 0.095 m^s"^, 
and 0.057 m^s^^ respectively. These values generally fall between the JONSWAP average 
value r = 0.038 m^s“^ , and the one obtained from observations in the Great Australian 
Bight (Young & Gorman, 1995), F = 0.152 m^s~^. Although the JONSWAP formulation 
with the widely used value F = 0.038 m^s~^ gives reasonable wave height predictions in 
most conditions (figures 13, 15 and 18), it cannot reproduce the observed variations in swell 
decay, and significantly overestimates wave heights in active ripple conditions, as was also 
noticed by Weber (1991a). In contrast the constant roughness = 4 cm) proposed by 
Weber (1991) yields values of that are still too high (by 30 cm) in the October 19 case, 
but too low (by 2.5 m) in the November 18 case (not shown). The movable bed model, 
adopted from Tolman (1994) without any adjustments, captures this variability, and fine 
tuning of all the empirical parameters, should further improve swell predictions. However, 
the movable bed parameterization requires site-specific sediment data that are not always 
available. Without such data, operational wave models may be better off with more robust 
dissipation models (e.g. Weber, 1991a; Tolman, 1994; Luo & Monbaliu, 1994; Young & 
Gorman, 1995). 

2. Model efficiency 

The new hybrid Eulerian-Lagrangian model CREST, presented here, was used to 
investigate the effects of a movable sandy sea bed on the transformation of swell across a 
continental shelf. Other physical processes such as wave generation, resonant non linear 
interactions between waves (see for example Herterich and Hasselmann, 1980) and res¬ 
onant Bragg scattering of waves by bottom features (chapter III) can be incorporated as 
additional source terms in the energy balance equation. Hence the present model provides 
an alternative to the Eulerian finite difference scheme models commonly used for wave 
prediction. With Aat a typical number of time steps for a ray bundle, and a typical num¬ 
ber of interpolation coefficients and B"; for a given time step, the CREST wave model 
requires memory space for storing the interpolation coefficients that is a factor Aaj x W (of 
the order of 200 in the calculations presented here) larger than the space used for storing 
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the spectrum. Thus CREST requires much more memory per grid point than an Eulerian 
model, that only needs to store the spectrum. The hybrid approach is attractive for appli¬ 
cations where the spatial scales Ls of variations in the source terms are much larger than 
spatial scales Lr of refraction effects. An Eulerian model describing refraction with a finite 
difference scheme in space requires a grid resolving both Lr and Lr whereas the Eulerian 
grid in the present model needs to resolve only Lr because Lr is resolved by the precom¬ 
puted rays. If Ls is much greater than Lr this property reduces drastically the number Agp 
of grid points required for an accurate integration of the energy balance equation. Reducing 
Agp has the added benefit that in coarser grids fewer grid points are used to interpolate of 
the source term for a given ray bundle, thus reducing the number Ns of interpolation coef¬ 
ficients per time step. The implementation of an Eulerian finite difference scheme with a 
resolution of about 500 m would have similar memory needs as the calculations presented 
here. 

The considerable memory burden imposed by the storage of the ray information 
can be reduced by dividing the model domain into subdomains. The use of subdomains is 
clearly a compromise between a pure Eagrangian advection scheme and practical consid¬ 
erations. At the internal boundaries it re-introduces numerical diffusion in the advection 
and re-couples the rays to the grid for x > 0. Although not necessary in the application 
presented here it seems unavoidable for implementations of CREST on very large areas 
( e.g. > 10^ km^). Eurther economy on the computer resources can also be achieved by 
specifying a maximum number of time steps along the rays after which the energy is inter¬ 
polated from the neighboring grid point. This scheme also increases numerical diffusion. 
Eor one-timestep-long rays the numerical scheme becomes a finite bandwidth version of 
the piece-wise ray methods. 

The representation by refraction alone of the effects of small scale bottom irregu¬ 
larities is cumbersome in the present model, and may not reflect the entire complexity of 
that process. A statistical representation of the interaction of waves with the smallest scale 
bathymetric features (e.g. the wave-bottom Bragg scattering source term described by Has- 
selmann (1966), and Eong (1973) appears attractive because it would improve the physical 
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description of wave-bottom interactions, and the rays computed on smoother bathymetry 
would be less scattered, thus requiring a smaller number Ns of interpolation coefficients. 
This improvement is described in chapter III. 

The small number of grid points in CREST is also advantageous for complex source 
terms (e.g. quartet interactions between wave components) that are prohibitively expensive 
to evaluate accurately on a high-resolution grid. Furthermore the flexible model grid gen¬ 
erated from any arbitrary set of points can be tailored to the bathymetry and shape of the 
model domain with higher resolution in the shallowest parts of the domain. In this respect 
CREST is similar to the TOMAWAC model (Benoit et al, 1996). 

For practical applications, computing ray trajectories is too expensive to allow a 
time-dependent ray geometry. This prevents the use of CREST in regions with strong 
temporal medium variations such as unsteady currents and tidal depth changes found in 
shallow estuaries and macrotidal seas, unless an approximate representation of these effects 
as source terms is found. 

G. SUMMARY 

A non-stationary spectral wave model was developed using a hybrid Eulerian-Ea- 
grangian scheme to examine the damping of swell propagating across a wide, shallow 
continental shelf. The model accurately represents refraction by advecting wave energy 
from deep water along a full spectrum of precomputed ray trajectories to a large number 
of grid points on the shelf. The source term in the energy balance is computed at each 
of these grid points, based on the complete frequency-directional spectrum. The source 
term is then interpolated from the Eulerian grid onto the rays, thus allowing for nonlinear 
coupling of wave components traveling along different rays. The energy balance is aver¬ 
aged over ensembles of rays to represent a finite spectral bandwidth. The (Eagrangian) 
computation of energy advection along rays and the (Eulerian) source term evaluation are 
carried out in parallel through the entire model domain. Source term formulations can be 
adapted from existing third-generation wave prediction models, whereas the finite differ¬ 
ence propagation schemes of these models are replaced with a full Eagrangian ray method. 
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This hybrid scheme avoids the numerical diffusion and ‘garden sprinkler’ problems of ex¬ 
isting models that use finite difference schemes. The ray calculations and source term 
interpolation scheme add considerable computational effort, but both the ray trajectories 
and interpolation coefficients are precomputed for a given coastal region and model grid. 
The spectral energy balance is integrated in time with an efficiency comparable to existing 
finite-difference schemes. 

The model was implemented with a source term restricted to energy dissipation in 
the bottom boundary layer over a movable sandy bed, as parameterized by Tolman (1994). 
The model was used to hindcast swell evolution across the North Carolina continental shelf 
for a range of wave conditions (significant wave heights between 0.5 and 10 m, and peak 
periods between 8 and 17 s) observed during two storms in 1994. Good agreement between 
observed and predicted variations of significant wave heights and mean wave directions 
across the shelf supports the hypothesis that refraction and movable bed bottom friction 
dominate the evolution of swell over a shallow sandy continental shelf. 
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III. BRAGG SCATTERING OF RANDOM WAVES 
BY BOTTOM IRREGULARITIES 

The energy balance (1.3) implicitly assumes that depth variations are small at scales 
comparable to the wavelength. Topographic features at those scales should therefore be 
excluded from refraction calculations but their effect on the waves can be represented in the 
energy balance in the form of Bragg scattering source terms. This stochastic representation 
of energy exchanges between wave components, was first developed by Hasselmann (1962) 
for wave-wave scattering, and later extended (Hasselmann, 1966) to wave scattering by 
external perturbations using a decomposition of the bottom elevation in Fourier modes. 
This theory is extended in the present chapter to heterogeneous wave and bottom properties, 
using decompositions in slowly modulated Fourier modes. Only the lowest order wave- 
bottom Bragg scattering process is considered. The present chapter reproduces with more 
details an article currently in press in the Journal of Fluid Mechanics (Ardhuin & Herbers, 
2001 ). 

A. SCATTERING THEORY EOR RANDOM WAVES IN HET¬ 
EROGENEOUS CONDITIONS 

The present derivation of the energy balance equation for random waves propagat¬ 
ing over an irregular sea floor uses a perturbation expansion of the wave energy, closely 
following Hasselmann’s (1962) derivation of energy transfers in quartet wave-wave inter¬ 
actions, and a ray approximation of medium variations adapted from Mei (1989, ch. 3). The 
result is a local energy balance equation that incorporates refraction and shoaling by large 
scale depth variations, and a source term describing Bragg scattering by seabed topography 
with small horizontal scales (of the order of the surface wavelength). 

I. General formulation 

We consider weakly nonlinear random waves propagating over an irregular bottom 
with a slowly varying mean depth and random small scale topography. For the sake of 
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Z = - 


Figiue 19. Random waves and iiiegular bottom: definition sketch. 


H(x) 


simplicity we will neglect the effects of mean ciments on wave propagation (see for ex¬ 
ample Bretherton & Gaiiett, 1969) and on wave scattermg by bottom imdulations (Kuby, 
1988; Ting, Lin & Kuo, 2000). All variables are non-dimensionalized with a representa¬ 
tive wavenimiber ^o, gravity acceleration g and water density p. Tire bottom elevation is 
represented by r = —H{x) +h{x), where h is a zero-mean small deviation fiom the gen¬ 
tly sloping large scale bottom features represented by —H{x), x is the horizontal position 
vector, and r is the elevation relative to the mean water level. The vertical position of the 
ocean fiee siuface is given by ^{x,t) with a zero mean value (figiue 19). Assiuning mo- 
tational flow for an incompressible fluid, the horizontal velocity field u is equal to V(t), the 
horizontal gr adient of a velocity potential, and the vertical velocity w is equal to 3(|)/dr. We 
fiulher assume that p is constant. Tlie goveniing equations for (|) are 

V2(})+^ = 0 for + (in.l) 


3(1) 

— = 'V(^-(yji— VH) at z = —H-\-Ji, 
oz 


C- 


(ni.2) 


(m3) 


3/2 ' 3r 3/ 3r3/3r 

(Hl.l) is Laplace’s equation, (III.2) is the ‘fiee slip’ bottom boimdary condition, and the 
‘combined’ smface boimdary condition (in.3) is obtained by eliminating linear teims in¬ 
volving ^ from the dynamic (i.e. Bernoulli’s equation) and kinematic conditions at the free 
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surface (see for example Hasselmann, 1962). ^ is given by Bernoulli’s equation, 


''"■•(I 


at z = ^. 


(III.4) 


Assuming that h varies on scales of the order of the surface wavelength, we in¬ 
troduce three small parameters: the wave slope £ = koUQ, the small scale bottom slope 
T] = koho, and a measure |3 of the large scale bottom slope |V//|. (IIL1)-(IIL3) are scaled 


V^(|)-1-^=0 for {-H+ r\h) <z<e^, 


(IIL5) 


^ = V(|) • (riVh - pv//) at z^-H + r\h, 
az 

3^(|) ^^, 3V(|) d(|) d^(|) 

^ + ^ = eV(^-V^-£V(^-—-at z = e^. 

dt^ dz dt dz dtdz 


(IIL6) 


(III.7) 


Following Keller (1958) we introduce slow space x = ax and time t = yt variables, h and 
(|) are assumed to be semi-stationary random processes in horizontal space and time (for (|) 
only), with evolution scales (aA:o) ^ and Y ^ respectively (Priestley, 1965), that can be 
decomposed into Fourier modes with slowly varying amplitudes. Following Hasselmann 
(1962) we shall approximate h and (|) with discrete sums, and take the limit to continuous 
integrals after deriving expressions for the evolution of the phase averaged wave energy. 
We write 

/^(x)=£Bl(^)ei'•^ (III.8) 

1 

where 1 are regularly spaced wavenumbers of bottom undulations and Bi are slowly varying 
amplitudes. Anticipating the effects of refraction, (|) is decomposed as 


k 


(III.9) 


where k are regularly spaced surface wavenumbers, and each k-component has an ampli¬ 
tude 4>ii, an eikonal Sk, and a local wavenumber 


k, (k, px) = V5k (x) 


(III. 10) 


such that k^ = k at the origin x = 0. 4>k and k^ are Lagrangian variables following a wave 
component along a ray trajectory. The spectral decomposition (III.9) for an evolutionary 
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process is ‘unique’, in a sense defined by Priestley (1981, theorem 11.2.3), only for a finite 
region in spaee and time, and is used here only to evaluate loeal variations of 4>ii. 

The slow spatial variations of 4>k ean result from shoaling, refraetion, and seatter- 
ing proeesses, as well as non-stationary and non-uniform wave eonditions. Sinee (|) and h 
are real it follows that k and Bi = B \, where the overbar denotes the eomplex 

eonjugate. 

In the vieinity of x = 0 the deeomposition (III.9) reduees to a Fourier sum 

= + (III.ll) 

k 

The simplified deeomposition (III.ll) will be used when no spaee differentiation is in¬ 
volved, taking advantage of the orthogonality of Fourier modes. 

The goal of the present derivation is to determine from (III.1)-(III.4) the energy 
balanee at x = 0 for eaeh k-eomponent of the wave speetrum (III.9). The solution depends 
on the relative magnitudes of the five small parameters: a, (3, y, p, and £. Here we use 

(III. 12) 

The ehoiee of a small seale bottom slope p mueh larger than the large seale slope (3 is 
usually well suited to sandy eontinental shelves, with the exeeption of the steeper beaeh 
and shelf break regions. This ehoiee makes the present derivation a priori different from 
Mei’s (1985) theory in whieh (3 ~ Tj. 

Following the method of Hasselmann (1962), the solution to (III.5)-(III.7) is ob¬ 
tained through a perturbation expansion in powers of £, 

(|) = (|)i-h £(|)2 + £^(|) 3 + h. o. t. (III. 13) 

The boundary eonditions (III.6) and (III.7) are expressed at z = —H and z = 0, respeetively, 
using Taylor series expansions of (|) about z = —H and z = 0, e.g. at the bottom, 

+ h. o. t. (III. 14) 

z=-H 

Eaeh term in (III. 13) will be found to be of the form 

^ (krZ + krH) giSi,k(x) bound wave terms , (III. 15) 

k s yKyti j 


Mz=-H+h — 


d(|) 

dz 


2 ^ 5 

7 =- h ^^ 2 az2 
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where kr is the magnitude of the local wavenumber vector k^, 5 is a sign index (+ or —), 
is the amplitude of the free wave component (k, 5 ) that propagates in the direction of 
5k„and^ = d>rlk. 

The slowly evolving spectral statistics of free wave components can be expressed 
in terms of the covariances F^- of the velocity potential amplitudes: 








(HI. 16) 


where the angular brackets denote an average over many realizations of the wave field, 
and in local space and time over a region that is large compared to the ‘fast’ scales of 

sea surface excursions but small compared to the slow scales (ako) ^ and (gko) of 

spectral variations. The contribution of the complex conjugate pairs of components (k, +) 
and (—k, —) are combined in (111.16) so that is the covariance of waves propagating 
in the direction of k. Note that the wavenumber separations ■k^= (Akr^xi ^kr,y) the sum 
(III.9) vary along rays owing to refraction. In the limit of small wavenumber separation a 
continuous cross-spectrum can be defined at x = 0 (e.g. Priestley, 1981 ch. 11) 


(0,t,k) = lim 
^ ^ iBkHo AkxAky 


(III. 17) 


The definitions of all spectral densities are chosen so that the integral over the entire 
wavenumber plane yields the total covariance of (|),- and 

The slowly varying bottom elevation spectrum in discrete form is given by Ff = 
{BiB i) and in continuous form by 


F®(x,l)= lim 

IrtHo AixAiy 


so that 


{h^ (x)> = 


-00 /»-[-00 


F«(i,l)dlxdl 


y 


(III. 18) 


(III. 19) 


This definition differs by a factor 2 from the one chosen by Hasselmann (1966) and Long 
(1973). 

The total wave energy at x = 0, in non-dimensional form. 


F(0,0 = 





1 


dz) + -{^ >, 


(III.20) 
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can be written as 


£ (0, t) = r r [e^E2 (k) + £3^3 (k) + £4^4 (k)] dk^clky + O (£3) , (IIL21) 

where 

£2(k)=£i,i(k), (IIL22) 

£3(k)=£2,i(k)+£i,2(k), (IIL23) 

Ea (k) = £ 2,2 (k) +£ 3,1 (k) +£ 1,3 (k) . (III.24) 

Here j (k) is the (/ + 7 )* order energy contribution from correlations between and 

yth Qj-(jgj- components with wavenumber k. Since the average in (111.20) is over several 
wavelengths, correlations between wave components with different wavenumbers that re¬ 
sult from reflections (i.e. standing wave patterns of nodes and antinodes) are averaged out 
and do not contribute to (111.21). For all (/, j) pairs Eij (k) = Ej^ (k). Hasselmann (1962) 
discarded odd-power energy terms £3 and E^ under the assumption that the sea surface is 
Gaussian. It was later found that this assumption is unnecessary (Benney & Saffman, 1966 
; Newell & Aucoin, 1971) as dispersion decorrelates the wave components during their 
propagation. Here additional terms involving correlations between two wave and one bot¬ 
tom component contribute to £ 3 , but these terms are shown to be bounded in appendix A. 
The dynamically important growing terms will be found in the 4* order energy E 4 (III.24). 

For freely propagating waves the potential and kinetic energy contributions to (III. 20) 
are equal and Eij is approximately given by the linear relation 

Eij (0, Ek)= kE^j (0, t, k) tanh {kH ). (III.25) 

Neglected in (III.25) are the contributions to the kinetic energy integral (III.20) from the z- 
intervals [—El -f- /z, —H] and [0, . Although these contributions are O (£^) and thus should 
be included in £ 4 , their magnitude is bounded and thus their time derivative is O (£^). 
All O (£^) bounded terms resulting from the surface and bottom boundary conditions can 
be discarded in the following analysis of energy transfers within the wave spectrum (see 
Hasselmann, 1962, for a detailed discussion). 
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2. First order solution 

Substitution of the first order wave field (|)i 


y, cosh(t,;+ <:,//) AgiSiiW 

cosh(t,//) 

in the surface boundary condition (III.7) yields 


(111.26) 


(111.27) 


where the radian frequency (0 {k) is constant along rays, and is given by the linear dispersion 
relation (in non-dimensional form): 


(Si{k) = [A:^tanh(A:^//)]2. 


(III.28) 


The slow space and time modulations of and the associated variations of the energy 
spectrum E 2 (k^) are not constrained by the first order equations, but can be determined 
from the fourth order energy E 4 (k^) (III.24), that depends on both second and third order 
waves. 

3. Second order solution 

Substituting (III. 13) in (III.5)-(III.7) and collecting terms of order £ and p yields 
the governing equations for the second order velocity potential (|)2 


V >2 + ?t^=0 for H<Z<0, 

az^ 

(III.29) 

a^2_ +V(^i-Vfi at z= H, 

az oz^ 

(III.30) 

+ at a = 0, 

oz 

(III.31) 


where NL 2 contains the non-linear terms in the surface boundary condition that force a 


bound wave solution (ItS* (Hasselmann, 1962, (47)). Note that refraction and shoaling terms 


associated with the large scale bottom slope V// are of higher order and do not contribute 
to the second order equations. Therefore ray curvature effects on (|)i can be neglected, and 
we can use k^ ^ k and 5k (x) ^ k ■ x in the vicinity of x = 0. A general solution to Laplace’s 
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equation (111.29) is formed with a Fourier sum of free and bound wave components with 
amplitudes 4>2and ^2 k' 


k,j 


cosh {kz + kH) 
cosh [kH) 


^lkW + 


sinh(fe + fc//)^si,s,,^ 
cosh{kH) 



(111.32) 


d> 2 follows from substituting the first order wave field (III.26), (III.27) in the right hand 
side of the bottom boundary condition (III. 30). 


^zkW = -E‘^:‘^«k k'q>i,k'e 




(III.33) 


k',j 


where (( 0 ^k') obey the dispersion relation (III.28). The bound wave 4>2k effectively cou¬ 
ples the bottom and surface waves. Substitution of (III.32) and (III.33) in (III.3I) yields a 
forced harmonic oscillator equation for the free wave amplitude 4>2 

A + 0)2^ {t) = E tanh {kH)] !i^5k_k'^pk' (0 • (ni.34) 

Following the method of Hasselmann (1962), the time derivative of the energy density 
El,! (k) of the second order waves in the limit of large t at x = 0, can be written in the form 
(appendix A) 

=i^(/t,//) ^^"cos 2(0-0')F® (k-k')£2 (k') d0', (III.35) 


where k = (A:cos0,A:sin0), k' = (A:cos0',A:sin0'), and 


K{k,H) 


sinh {2kH) [2kH -f sinh {2kH)] 


4. Third order solution 


(III.36) 


Slow modulations of (|)i yield third order terms in Laplace’s equation. Substituting 
(III.13), (III.26), and (III.27) in (III.5)-(III.7), collecting terms of order E^, er\, ri2,a, |3, 
and y, and using the approximations (in the vicinity of x = 0) k^ = k + O (|3x), and 5k (x) = 
k ■ X -|- O (ax, (3x), yields the following equations for the third order velocity potential (|)3 

I 


V2(^3 + 


dz^ 


-i y k ■ V f $1 ^Q^h(fc,Z + fc//) ^ ^i(k.x-sa)t) 
^ V ’ cosh{krH) 
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II 


-iVV-fk (^rZ + krH) \ i(k.x-sa)t) 

ts V ' cosh (krH) ) 

for -H<z<0, (111.37) 


dz 



V 


^ + V (^2 ■ V/i -i £k ■ 

k,i 


i(k-x—scot) 


at z = — 


H, 


(111.38) 


VI 


^ ^^i^^^(^^^ei(kx-scot)^NL 3 at z = 0, (III.39) 

dt^ dz dt 

Note that third order terms involving (|)i in the bottom boundary condition (III.38) vanish 
because d^(|)i/dz^ = 0 and d^i/dz = 0 at z = —H. The right hand side forcing terms of 
(III.37-III.39) include Bragg scattering terms (III and IV), effects of spatial heterogeneities 
(I,II, and V), non-stationarity (VI), and third order nonlinear surface terms that are gathered 
here in the term VL 3 . This set of equations is linear in (|) 3 . Therefore (|)3 is the sum of a 
homogeneous solution (absorbed in (|)i) and four particular solutions. 


4*3 = 4'f + 4'3‘' + 4'^+4'3^ 


(III.40) 


where sc, he, ns and nl, stand for scattering, heterogeneity, non-stationarity, and non¬ 
linearity, respectively. Each solution satisfies (III.37)-(III.39) forced respectively by the 
scattering terms (III and IV) only, the heterogeneity terms (I, II and V) only, the non- 
stationarity term (VI) only, and the surface non-linearity terms (VL 3 ) only. Although (j)"* is 
resonantly forced, it contributes only bounded terms to E 4 (Hasselmann, 1962). Similarly, 
nonlinear contributions to the scattering terms III and IV (the O (e^T]) products involving 
and bottom undulations) yield only bounded contributions in E 4 . The remaining solu¬ 
tions (1)1^^, (|) 3 ®, and i])"® contribute growing terms, , and to E 4 . Following the 

method used to obtain dE^' /dt, at x = 0 we have (appendix A) 


£|'l(k)+£f^3(k) 


dt 


fin 

K{k,H)J^ cos2(0-0')E®(k-k')E2(k) d0', (III.41) 


d 


E3 -(k)+E-(k) 


dt 


dE2ik) 

dt ’ 


(III.42) 
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£3‘'i(k)+£5(k) 


he 


dt 


= -c,(k)-v£2(M, 


(IIL43) 


where Cg (k) is the group veloeity of linear waves (A.27), and E 2 (k^-) is a Lagrangian 
variable that deseribes energy evolution along the ray trajeetory [x(k,(3r),k^(k,Pr)] of 
wave eomponent k, where r is the along-ray eoordinate. E 2 (k^) is defined by 


Fij(x,t,k;.) = lim 

lakrl-^O Akr^x^^r,y 

E 2 (x, t, k^) = KE^i (x,t, k^) tanh {krH). 


(111.44) 

(111.45) 


Note that the adveetion term Cg (k) ■ VF 2 (k^) deseribes the divergenee of the energy flux 
in Lagrangian eoordinates, and thus incorporates refraction and shoaling effects. All other 
terms in (III.42)-(III.43) depend only on the energy at x = 0 where the Lagrangian wavenum¬ 
ber kr is equal to the Eulerian wavenumber k, and E 2 (k^) = E 2 (k). 


5. Energy balance 


Combining (III.35), and (111.41 )-(III.43), the rate of change of the fourth order 
spectrum (111.24) at x = 0 is given by 


3^4 (k) 

~^r~ 


-^^-Cs(k).V£2(M 

r2n 

+K {k, H) 1^ cos2 (0 - 0') F® (k - k') [F 2 (k') - E 2 (k)] d0', 

(III.46) 


where K{k^H) is given by (111.36). 

To assure that F 4 is bounded for large t, that is dE^/dt = O (e^) , the right-hand side 
terms of (111.46) must balance. Recognizing the first two of these terms as the total deriva¬ 
tive of E 2 (k^) along a ray trajectory, and replacing E 2 by F, we obtain (using dimensional 
and unsealed variables from now on) the Lagrangian energy balance equation at x = 0, 
accurate to order £^: 

= ^Bragg (k) + O (e^) , (III.47) 


'■271 


5Bragg(k) =47tg2//^2%(/t//) / cos^ (0 - 0') F® (k - k') [F (k') - F (k)] d0', (III.48) 
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Figiue 20. Values of x {kH), as defined by (111.49). 


with 




{kH)^ [tanh(A'/f)]5 


(m.49) 


siuh {2kH) [2kH + sinh {2kH)] 

(III.47) describes the net energy transfer at x = 0 to a wave component with wavenumber 
k (propagating in diiection 0), resulting fiom triad interactions involving a wave of the 
same radian fiequency (O and a different wavenumber k' (diiection 0'), and a bottom com¬ 
ponent with the difference wavemmiber 1 = k — k' (figiue 3). The energy transfer between 
components k and k' is proportional to the energy difference of the wave components and 
the bottom spectmm density at 1 = k — k'. The factor cos^ (0 — 0') in (III.48) indicates that 
there is no energy hansfer between waves propagating in peipendicular diiections. The fac¬ 
tor X {kH) has a single maximum, approximately equal to 0.049 for the mteiinediate water 
depth kH ^ 1.27 (figiue 20). In addition to duectional and waveniunber dependencies, 

_9 . 

the scatteiing sfiength is proportional to .^7 ^, mcreasing strongly with decreasing water 
depth. Taking into account then different noimalization of the bottom elevation spectmm. 
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the present expression (111.48) of Ssragg is 4 and 8 times smaller than the expressions given 
by Long (1973) and Hasselmann (1966), respectively. 

6. Conditions of validity 

The present theory is both a spectral generalization and higher-order energy con¬ 
serving form of the solution given by Davies (1979) for sinusoidal bed undulations. Davies 
describes the generation of second order waves (|) 2 , but uses constant amplitudes for (|)i, and 
thus does not account for the associated energy losses of the primary waves. In the present 
theory the extension of the perturbation expansion to third order provides the balancing 
terms Eys and £'31 (III.41), necessary for the conservation of the total energy in (III.47). 
Whereas Davies’ theory assumes small reflected wave amplitudes, (III.47) can describe 
finite cumulative reflections over large distances and even complete localization of waves 
over rough bottom topography. However the present theory assumes that significant wave 
amplitude variations occur over scales of 0 {a) wavelengths with a ~ £^, and thus cannot 
accurately describe strong localized scattering that modify the wave amplitudes over scales 
of only a few wavelengths (see Mei, 1985, for a discussion of those effects over sinusoidal 
bottom topography, including in particular the importance of near-resonant interactions in 
that case). 

It should be noted that the wavenumber spectrum E {k) = jQ^kE (k)d0 was as¬ 
sumed continuous in order to derive (III.47) in the limit of large times, removing the sin¬ 
gularities for perfect resonance in (A. 6 ), and reducing the bandwidth of important near¬ 
resonant interactions to a region of the spectrum where E (k) can be considered constant. 
Thus (III.47) is not valid for monochromatic waves. Whereas the initial growth of the scat¬ 
tered energy is proportional to for resonant monochromatic waves, it is only proportional 
to t for waves with a continuous spectrum, because resonance becomes more selective with 
time, affecting a wavenumber bandwidth that narrows proportionally to 

Another consequence of the asymptotic large time limit taken in appendices A-D, 
is that the stochastic model (III.47) may not describe accurately wave evolution over natural 
sea beds which are often not homogeneous over scales of many wavelengths. The robust- 
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ness of (111.47) for short propagation distances is examined in § 3 through comparisons with 
deterministic models for wave evolution over a finite patch of sinusoidal bars. (III.47) in¬ 
cludes wave-bottom interactions with |k — k'| <^k, violating our scaling assumption I 
This particular aspect is discussed in § 3.4. 

7. Extensions of the present theory 

The present energy balance (III.47) may be extended to higher orders of t) and/or £ 
by closing the energy Taylor expansion at E(^, giving an evolution equation for £2 + £^ 4 - In 
the case of steeper waves, say a ~ |3 ~ Y~ ^ £^, it can be seen that all the energy transfer 
terms derived here (III.35, III.41-III.43) are moved from £4 to joining the additional 
source term 5ni that represents resonant quartet wave-wave interaction (Hasselmann, 1962; 
Herterich and Hasselmann, 1980). Extensions to steeper waves and steeper topography, 
for example a ~ |3 ~ y ~ ^ £^, should yield at least two additional source terms, cor¬ 

responding to higher order Bragg scattering (class II and III, see for example Liu & Yue, 
1998). 

Furthermore it can be expected that including higher order heterogeneity effects 
and non-linearity should introduce nonlinear effects on the left hand side of (III.47), as 
described by Willebrand (1975). For example in the present theory E^ contains correlations 
between the tertiary waves t])"' and the heterogeneity and non-stationarity terms (j)^® and t])"® 
. Thus it may be possible to derive a more complete energy balance equation with not only 
the source terms for the individual physical processes that contribute to the evolution of the 
wave spectrum, but also the cross-interactions of these processes that are usually neglected 
in wave prediction models (Komen et al, 1994). 

B. RANDOM WAVES OVER SINUSOIDAL BARS 

Following Davies (1979) we consider a simple seabed consisting of sinusoidal bars 
on an otherwise flat bottom for which analytical results exist that have been verified in 
laboratory experiments. Waves arriving from x — — 0 ° at an incidence angle 0/ are partially 
reflected, in a direction 0 i; = 7 t — 0 /bya patch of m sinusoidal bars of amplitude b, aligned 
with the y axis. The barred profile h = b sin (lx) covers the region —L < x < L where 
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Figure 21. Wave scattering by sinusoidal bars 
Schematic of incident waves (dashed crests) and reflected waves (dotted crests) on a patch 
of sinusoidal bars (gray shades). 

L — tniz/l and / is the bar wavenimiber (figiue 21). The incident wave field is assumed 
to be a coutmuous spectnun Ej (k) of uniduectioual (0 = 0/) waves. The total reflected 
energy Er in the far field (.y <§; —Z) predicted by the stochastic and detenninistic theories 
are compared for both normal and oblique incidence cases, in the limit of large ni and for 
finite Jii. 

1. Stochastic source term approach 

In a steady state imifoim along the y-axis, the energy balance (III.47) for the bottom 
profile described above simplifies to 

dE(k) 

CgCOS0— -—- = ^Bragg (k) . (111.50) 
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Figme 22. Resonant wavenumbers for sinusoidal bars. 

Resonant triads for waves over sinusoidal bars represented on the wavenumber plane. The 
bar wavenmnbers are fixed at (/,0) and (—/,0) and all possible pairs of resonant smface 
waveniunbers ki and Icr lie on the vertical dashed lines. For the wave duection 0 shown 
here, the directional spectral density E (0) (the integral of kE (k) along the thick anow) is 
affected by energy transfers m the resonant (k/, k^, 1) triad. 


In order to evahrate Ssragg (111.48) we approximate the finite patch of sinitsoidal bars as 
a subsection of a sinirsoidal bottom extending to infini ty, for which the bottom variance 
spectnrm is a doirble Duac distribrrtiorr 

^^0) = ^ [5(/,0)+5(-/,0)] . (m.51) 


Orrtside the barred section (|.r| > L) F^ is set eqiral to zero. The singirlarity in (III.51) is 
removed in ^Bragg by mtegrating (in.50) over k for a fixed duection 0 (Figrue 22). Changing 
variables from (A', 0') to the corresponding resonant bottom wavenirmber 
(Ixjy) = A'(cos0 —cos0',sin0 —sin0') we obtain 


dE(0) 

dx 


where 


J]L f f cos^ (9 - (1) [£ (k - 1) - £ (k)| ^ 

cosey 7k.i>0 [l_cos(e-e'))[2i-//+smh(2«/)]^ ' 

roo 

E{Q)= / A£'(A'cos0,A'sin0)dk 


(m.52) 


(m.53) 
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is the directional spectrum integrated over all wavenumbers, and the Jacobian 
J = 1/{k[l — cos(0 — 0')]} of the transform from (A:,0') to {Ix^ly) is used. Note that the 
integration over 1 is restricted to the half plane where k ■ 1 > 0 (figure 22). 

For —n/2 < 0 < 7t/2 (III.52) describes the evolution of an incident component 
with direction 0/ = 0. Only interactions in the neighborhood of the resonant triad k/ = 
/(l,tan0)/2, k^ = / (—l,tan0) /2, 1 = (/,0) contribute to this integral (figure 22). For 
7t/2 < 0 < 37t/2, (III.52) describes the evolution of a reflected component with direction 
0^ = 0 resulting from the resonance of k/ = / (1, — tan0) /2, k/; = / (—1, — tan0) /2, 1 = 
(—/,0). Substitution of (III.51) in (III.52) yields 

^ = -D;,[£(k/)-£(k^)] for -L<x<L, (III.54) 


where 


Dx = 


7r(?^cos^ (20/) 


J^ + sinhl 

cos 9 / T^^inn I ^.05 0^ 


(III.55) 


16cos^0/ [1 +COS (20/)] 

For weak reflection (E (k/) 3> £ (k//)) we can neglect changes in £ (k/). Integrating 
(III.54) from L to —L yields 


Er = ILDxE (k/) for a: < —L. 


(III.56) 


For unidirectional incident waves with a spectrum £'/(k) = 5(0 —0/)£'/(k) /k, the total 
reflected energy is given by 

Er = DIE, {k ,), (III.57) 


where D is an non-dimensional coefficient 

2LDx rnn^b^ cos^ (20/) 


D = 


Ik, 


4cos^ 0/ [1 + cos (20/)] 


IH 

COS0/ 


-l-sinh 


IH 

cos 6 / 


For the particular case of normal incidence (0/ = 0) D reduces to 

D= - 

S [IH + smh {lH)f 


(III.58) 


(III.59) 
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2. Comparison with deterministic theory for normal incidence 

In Davies’ (1979) theory for weak refleetion of a normally incident monochromatic 
wave train by a patch of m sinusoidal bars with amplitude b, the ratio of the reflected and 
incident wave amplitudes is given by 

V 2M (-ir2tsm(2tt) 

™ 2*H + smh(2*H) / {ff-l 

Theoretical values of Kqh have been verified experimentally by Heathershaw (1982; see 
also Davies & Heathershaw, 1984), even in cases with large reflection coefficients. 

For random waves with a wavenumber spectrum Ei {k) the reflected energy £’^,dh 
is the convolution of | Kdh {k) \ and Ei{k), 


= / |kdh(^)| Ei{k)dk. 


JO 


(111.61) 


2 

The response function |kdh| has a ‘resonant lobe’ of width n/L and height proportional 
to centered at the resonant wavenumber k = l/2, and narrower side lobes at higher and 
lower wavenumbers (figure 23,dotted curve). In the limit of large m (equivalent to the large 


t limit in the stochastic theory), iKonf approaches a Dirac distribution 


I Kdh I 


m'Krl 


2 bk 


8 \2kH + smh{2kH)) 


Y/2k 


^ 15(2 


(III.62) 


Since Ej {k) is continuous, the substitution of (III.62) in (III.61), yields 




rnn^b^P 


S[lH + smh{lH)] 


,El[^ 


(III.63) 


which is identical to the stochastic theory prediction (III.57),(III.59). The exact agreement 
of the stochastic and (experimentally verified) deterministic theories in the limit of large m, 
where both are valid, confirms that the coupling factor % (III.49), which differs by factors 
of 8 and 4 from previous publications (Hasselmann, 1966; Long 1973), is correct. 


3. Oblique incidence and finite numbers of bars 

Davies’ (1979) theory for wave reflection from sinusoidal bars was generalized to 
oblique incidence and finite reflection coefficients by Mei (1985), using an approximation 
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kll 


Figure 23. Reflection of waves by sinusoidal bars 

'y 

Response fiinction IkdkI H = 25 In/l = 300 ni, 6 = 0.05 m, in = 4, and three 

incidence directions: ., 0/ = 0;-, 0/ = 60°;-, 0/ = 75.5°. (coiiesponding to 

resonant wavenumbers k = //2, k = l, and k = 21). The response fimctions are normalized 
by then maximum vahtes indicated on the figine. For reference a generic wave spectrum 
is inchrded (o, arbitrary units) with a Pierson-Moskowitz shape, a peak period of 14 s, and 
imiform infragravity energy levels. The total reflected energy is the convohrtion of |kdk| 
and the wave spectrum. 
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for weak detuning from resonance. Dalrymple & Kirby (1986) applied Mei’s theory to a 
finite patch of bars and derived the amplitude reflection coefficient Kdk (their equations 5 
and 9). For normal incidence Kdk is in good agreement with the experimental results of 
Davies & Heathershaw (1984), and reduces to Kqh in the limit of small bar amplitude b. 
For oblique incidence no experimental verification exists but Mei’s theory was verified nu¬ 
merically with solutions of Kirby’s generalized mild slope equations (Kirby, 1993). Values 
of I Kdk I for a patch of four bars of wavelength In/l = 300 m, and amplitude b = 0.05 m 
in 25 m depth, are shown in figure 23 as a function oik/I for different incidence an¬ 
gles 0/. The interaction between the bottom undulations and the surface gravity waves is 
dominated by near-resonant triads, for which |kdk| is maximum. The Bragg resonance 
condition k = 2l/ cosO/ determines the wavenumber for which reflection is maximum, as 
a function of the incidence angle. For example, for the wave spectrum shown in figure 23, 
back scattering (0/ Ri 0, 0^; 180°) is confined to the low wavelength (infragravity) part 

of the spectrum, and shorter swells are scattered forward (see the response functions for 
0/ = 60°, 0^; = 120° and 0/ = 75.5°, Qr = 104.5° in figure 23). 

To determine the accuracy of the stochastic theory for a finite patch of bars, the 
total reflected energy Er predicted by (III.57)-(III.58), valid only in the limit of large m, is 
compared to the ‘exact’ F'/j^dk predicted by the deterministic theory (III.61 where Kdh is 
replaced by Kdk)> valid for arbitrary m. In these calculations Ej (k) is taken to be a Pierson- 
Moskowitz (1964) spectrum with a peak period Tp = 14 s. A white background spectrum 
E ik) = 0.04F'/ {kp) is added to represent contributions of longer wavelength infragravity 
waves (figure 23). The convolution integral (III.61) is computed numerically over the range 
0.005 <k/l < 4, for incidence angles 0/ = O, 60° and75.5°. Other parameters are 7/= 25 m 
and b = 0.05 m. A small b value was chosen to have a small reflection coefficient (kdk < 
0.1) because (III.57) neglects variations in the incident energy Ej and thus is valid only for 
weak reflections. The relative difference between stochastic and deterministic theories is 
shown in figure 24 as a function of m. This difference is sensitive to the variations of the 
wave spectrum across the resonant lobe and the relative magnitude of the side-lobes of the 
response function |kdk| • It vanishes in the limit m —> oo as the width of the resonant lobe 
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Figiue 24. Comparison of stochastic and deterministic theories 
Relative differences between the reflected energy predicted by the stochastic theory and the 
spectral form of the deterministic theory, (in.61, replacing Kdh by Kdk), as a fimction of 
the number of bars tu. All other parameters are the same as those used in figiue 23. The 
incident wave spectnim is shown in figiue 23. 

and the height of the side-lobes go to zero. As m mcreases the predictions of both theories 
converge, as expected since both theories are valid for large m. For all three incidence 
angles 0/ = 0, 60° and 75.5°, the difference in reflected wave energy predicted by the 
stochastic and deterministic theories is less than 25% for more than three bars. Tliis rapid 
convergence not only provides a fiuther consistency check of the couplmg factor x (in.49) 
for cases of oblique incidence, but also indicates that the stochastic Bragg scattering theory 
is siuprisingly robust. Although formally valid only in the asjmptotic limit of many bottom 
wavelengths, it yields reasonable estimates of energy transfers resultmg fiom scattering by 
only a few bottom imdulations. 
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4. Bottom slope effects 


The large time limit used to evaluate fourth order energy terms in appendices A- 
D requires implicitly that the large scale bottom slope does not significantly change the 
interaction over a distance Ar that allows waves to propagate across a sufficiently large 
number of bars to approach the asymptotic limit of the energy transfer (figure 24). 
Wave refraction by the large-scale bottom slope changes the surface wavenumbers and 
thus introduces a detuning of near-resonant wave-bottom interactions. This detuning effect 
can be neglected only if changes in the surface wavenumbers are small compared to the 
width of the resonant lobe of the response function |kdk| (figure 23). 

For simplicity we consider a finite patch of sinusoidal bars aligned with the y- 
axis, with wavenumber I, superimposed on a plane bottom with a downward slope |3 in 
a direction 0^. The along-ray gradient of the resonance mismatch u = (2A:cos0 —/) //is 
given by 


du 

dr 


2 

7 


7 Jk 


— A:sin0 



(III.64) 


Using SneTs law we have 


du 


-4|3A:^ cos 0ft 


(III.65) 


dr I [2kH -|- sinh (2kH)] 

For small bottom slopes the distance traveled by the waves across the bar field is Ar 
IniaTlI (Icosd), giving a change in the resonance mismatch 

27rmfl|3cos0ft 


Am 


cos^ 0 l 2 kH -I- sinh ( 2 kH)j' 


(III.66) 


Detuning of resonant interactions by refraction can be neglected if |Am| is small compared 
with the (normalized) width of the resonant lobe l/nta, that is 


nta |Am| 1. 


(III.67) 


(III.67) also follows from considering the phase difference between waves reflected by 
the first and m* bars, which should be small compared to k/2 to allow the constructive 
interference that causes resonance. 

(III.67) is a necessary condition for the application of the stochastic theory. For 
a given bottom slope |3, (III.67) imposes a maximum incidence angle 0max- For practical 
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purposes we assume that the largest acceptable value of |Am| is about 0.5, giving 


3 4nml^cosQb 

““ ^ 2kH + smh{2kH)' 


( 111 . 68 ) 


For example, considering 14 s period waves in 25 m depth with a bottom slope |3 = 2 x 10^^ 
at an angle db = 60°, and taking nia = 2, the source term (III.48) is expected to significantly 
overestimate the energy of scattered waves for incidence angles greater than Omax ~ 83°, 
corresponding to a ratio k/l = 4.3. 

It should be noted that (III.67) is consistent with the scaling of (III.1)-(III.4) requir¬ 
ing that bottom and surface elevations have comparable horizontal scales. This scaling is 
violated for large angle interactions (i.e. k 3> / for 0 close to 90°), even on a flat bottom. 
In the following the contribution of wave-bottom interactions to (III.48) is taken to be ac¬ 
curate for 0 < Oniax and is neglected for 0 > Omax- This crude truncation of the interactions 
is expected to give only qualitative results for the scattering of waves at large incidence 
angles. Sensitivity of predicted spectral evolution to the choice of the cut-off angle Omax is 
examined for natural shelf topography in § 4. 

C. HINDCAST OF WAVE SCATTERING ON A NATURAL SHELL 


The effect of Bragg scattering on directional wave spectra evolution is illustrated 
here with a numerical model hindcast of swell evolution observed across the North Car¬ 
olina shelf. The scattering source term ^Bragg (III.48) was implemented in the spectral 
model CREST (chapter II), that integrates the energy balance (III.47) in time using a hy¬ 
brid Eulerian-Lagrangian numerical scheme. In addition to ^Bragg a bottom friction source 
term Smc is included in the energy balance to account for energy dissipation in the boundary 
layer over a sandy movable bottom. Details of the model formulation, numerical scheme, 
treatment of boundary conditions, and parameterization of bottom friction are given in 
chapter II. ^Bragg is evaluated using bottom elevation spectra that were estimated from 
high-resolution bathymetry surveys. Processes not represented in the model such as wind 
generation, effects of currents, wave breaking, and non-linear effects, are expected to be 
negligible because at the time of the hindcast (21:00 Greenwich Mean Time, on October 
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20, 1994) local wind speeds (3 m s^^), and current velocities (~ 20 cm s"\ Lentz et al, 
1999) were weak, and the observed waves were long-period (~ 13 s) swell with low sig¬ 
nificant height 1 m). 

1. Wave data 

Frequency-directional wave spectra were estimated from measurements on the outer 
and inner shelf near Duck, North Carolina (figure 25 a). An array of pressure sensors, lo¬ 
cated 1 km from the shoreline in 8 m depth, was operated by the Army Corps of Engineers 
Field Research Facility (FRF), in Duck, North Carolina, and a 3 m discus pitch and roll 
buoy located close to the shelf break, in 49 m depth, was operated by the National Data 
Buoy Center (NDBC). Standard techniques (Berbers et al, 1999; § II.D) were used to ob¬ 
tain estimates of the frequency-directional wave spectra at both locations. The NDBC buoy 
wave spectrum was transformed across the shelf break to deep water using Snel’s law, as¬ 
suming parallel depth contours, in order to obtain the offshore boundary condition for the 
model. 

2. Bottom topography 

Bathymetric data for most of the shelf was available from the National Ocean Ser¬ 
vice (NOS). In regions not covered by the NOS archives, water depths were measured 
during instrument deployment and recovery cruises in a series of experiments (DUCK94; 
Sandy Duck; SHOWEX) on the North Carolina continental shelf, using a single precision 
depth recorder. Additionally, high-resolution multibeam sonar bathymetric surveys were 
conducted during the SHOWEX experiment in November and December 1999, in two 6 
km X 6 km square regions of the inner shelf (labeled 5i and S 2 in figure 25a). This data set 
was processed with the MB-System software (Caress & Chayes, 1995) to obtain 10 m res¬ 
olution grids shown in figure 25 {b,c). The vessel motion and tide were carefully removed, 
although a slight but systematic measurement bias is still noticeable in the striped pattern of 
figure 25 (c), yielding an artificial ridge of spectral densities on the x-axis of figure 26 {b). 
Although the high-resolution bathymetry data was acquired five years after the wave data, 
comparisons with depth soundings, performed within a few months of the wave data col- 
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Figiue 25. High resolution (10 m) bottom topography 
Bottom topography of the North Carolina continental shelf (a). The sqitares marked SI 
and S2 are the regions enlarged m (b, c). Other symbols indicate locations of the FRF 8 rn 
depth array (8M) and NDBC 3 m discirs buoy (44014). 
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Figme 26. Bottom elevation spectra 

Coutoiu' plots of the bottom variance spectra estimates for regions SI (o) and S2 (b). The 
contoiu’ values are logjo [4k~F^) with in m^rad“^, and the contoin mteival is 0.5. 
Circles indicate the bottom components that interact with waves arriving from the east with 
freqrtencies 0.05 (umer circle), 0.12 (rrriddle cfrcle) and 0.25 Hz (outer circle). Axes imits 
are reciprocal wavelengths Ix/ {2k) and ly/ {2k). (c) duection-mtegrated spectra for SI (— 

—) and S2 (-). The vertical lines indicate the bottom scales responsible for scattering 

0.08 Hz swell, for various incidence angles 0/ . 
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lection, show good agreement, suggesting that bottom topographic features of scales larger 
than 500 m have not moved in regions Si and S 2 . 

Bottom elevation spectra F® ( 1 ) and Ff ( 1 ) (figure 26a,b) were estimated for regions 
Si and S 2 , respectively, based on bidimensional Fourier transforms of Hanning windowed 
1.6 km X 1.6 km square segments with 50 % overlap. The large scale bottom slope was 
previously removed from each segment using a bilinear fit. Spectra of the large scale shelf 
topography (not shown), computed from the entire bathymetry grid, are consistent with 
the spectral levels at small I shown in figure 26a,b. The bottom elevation spectra are not 
isotropic, showing a preferential north-east/south-west orientation of intermediate scales 
features (200-1000 m) that are most important for swell scattering. It also appears that 
bottom spectral levels at these scales are about a factor 4 higher in region 5i (15-25 m 
water depth, variance 1.4 m^) than in the deeper region S 2 (20-40 m, variance 0.35 m^, see 
figure 26c). Lacking detailed topographic information in other regions, the bottom eleva¬ 
tion spectrum used in model hindcasts is taken to be uniform over the entire continental 
shelf. Hindcasts are presented in § 4 based on both estimates Ff ( 1 ) and ( 1 ), illustrating 
the likely range of scattering effects. 

3. Numerical model 

The numerical wave model CREST used for the present calculations is described 
in chapter II with unchanged spatial and wavenumber grids. The model consists of a pre- 
computation of wave rays and a Lagrangian time integration scheme for the energy balance 
(III.47). In contrast to more widely used finite-difference schemes (see for example the 
WAMDI group, 1988; Booij et al, 1999) the Lagrangian approach avoids numerical diffu¬ 
sion that could cause an artificial broadening of the wave spectrum in shallow water (not 
related to physical scattering processes). The Eulerian model grid, shown in figure 27 is 
unstructured and much coarser than the bathymetry grid. It consists of 329 points x„ dis¬ 
tributed over a large portion of the shelf between latitudes 35 and 37° N. Wave rays are 
traced backwards from the Eulerian grid points to the model boundary, using a smoothed 
(2 km scale) bathymetry grid that resolves wave refraction over the large scale shelf topog- 
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Figme 27. Model grid 

The giid points where the soiuce tem is evaluated are the nodes of the triangidar mesh. 
A linear interpolation is applied in each triangle to approximate the soiuce term along the 
rays. The 100 m depth contoiu is indicated by the dotted line. 
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raphy. Along each ray the energy balance (111.47) is integrated in time. At the grid points 
the full wave spectrum E (k) is evaluated using ensemble averages of rays within finite 
wavenumber bands / corresponding to 19 frequency bands fq spaced exponentially with 
a 5% increment from 0.05 Hz to 0.12 Hz, and 120 direction bands 0,- spaced linearly over a 
full circle with a 3 degree resolution. The spectral source terms Smc (k), representing bot¬ 
tom friction (§ II.B), and Ssragg (k), given by (III.48), are evaluated at the grid points based 
on the local spectrum E (k) and other parameters. Sfnc (k) and ^Bragg (k) are interpolated 
onto the ray trajectories to account for the energy losses (bottom friction), and exchanges 
with other wave components (scattering), of component k during propagation. Details of 
the time integration and interpolation schemes can be found in § II.3. 

The model was run here with uniform and steady offshore boundary conditions and 
a fixed integration time step At = 10 minutes, until a steady state was reached. To determine 
accurately the contribution of •S'eragg over the time step At, an implicit integration scheme 
was used. Omitting other source terms and propagation effects, (III.47) can be written in 
discretized form and for a given wavenumber magnitude k, as a set of linear equations 

= Angm~hikH) ik)E {k,dj) for all /, (III.69) 

where O,- are the discretized directions with k/ = A: (cos O,-, sin0/) and the matrix L(A:) is 
given by 


Lij ik) 


cos^ (Qi-Qj) E^ (k—k;) J^cos^(0,-e„)F^ (k, 

n 



A0, 


(III.70) 


with 5,y = 1 for i = jand 0 otherwise. As discussed below in § III.4, E^ (1) is replaced by 
zero in (III.70) for k/l greater than [k/l)^^-^. Since L is real and symmetric, it can be di¬ 
agonalized and represented as the matrix product L = VDV^ where D is a diagonal matrix 
with the eigenvalues X,- as diagonal elements, the columns of V are the corresponding nor¬ 
malized eigenvectors, and is the transpose of V. Using this decomposition the solution 
of (III.69) can be given in the form 


£(fc,e,-,t-hAt) = ^£l^-j(fc)exp A%gm-hikH)Xjik)^t Vij{k)E ik,Qi,t). (III.7I) 


j I 
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The source term Ssragg, R in (H. 15), is given by the average change in E {k, G,-, t) over a time 
step At 

^Bragg (k, 0,-) = [E {k, 0,-, t + At)-E {k, 0,-, 0] /At. (III.72) 

The matrices V {k) and eigenvalues Xi {k) are precomputed using Jacobi’s algorithm (see for 
example Press et al, 1992) for 500 values of k covering the entire range of wavenumbers in 
the model, and the resulting matrices V and D are interpolated on the spectral model grid. 
The high accuracy of the implicit numerical scheme was confirmed through comparisons 
with an explicit fifth order Cash-Karp Runge-Kutta scheme. 

4. Hindcast 

The model hindcast was performed both with and without the Bragg scattering 
source term, to isolate the scattering effects from other processes (refraction, shoaling 
and bottom friction), and using two different measured bottom elevation spectra (Ef (1) 
and ^^(1), figure 26) to estimate the possible variability of the scattering effects. Bot¬ 
tom components with wavelengths larger than 5 times the surface wavelength (k/l > 5) 
are excluded in the evaluation of Ssragg because, as discussed in § 3.4, the theory is not 
expected to be accurate for near-grazing angle interactions. Figure 28 shows an example 
wave spectrum predicted in region 5i (20 m depth, figure 25), and the corresponding Bragg 
scattering source term. The source term has a 3-lobe shape with negative values near the 
peak 0p of the directional wave spectrum, and positive maxima on both sides of the peak, 
at about 0p ± 30°. The interactions broaden the peak of the directional wave spectrum 
(forward-scattering) and cause weak, almost isotropic back-scattering. Sign reversals of 
SBragg within the main lobe (figure 2Sb) are caused by irregularities in the wave spectrum 
(figure 28a). Bragg scattering tends to smooth the directional wave spectrum, with an evo¬ 
lution time-scale E/S^ragg of the order of 10^ and 10^ s in 20 and 50 m depth, respectively. 

The combined effect of Bragg scattering and refraction is shown in figure 29 with 
the predicted cross-shore evolution of the mean wave direction (from) at the peak fre¬ 
quency, 0, taken as the direction of the first-order moment vector 

{ai,bi)= r(cos0,sin0) £(/p,0)d0, (III.73) 

Jo 
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Figme 28. Wave spectnim and Bragg scattering somce teim 
{a ) Predicted wave spectnuu E in region Si (figme 26) and (b ) conesponding Bragg 
scattering somce term Ssragg, based on bottom spectnim Ef. Contours in (Z> ) are solid 
for positive values (yellow to red color shades), heavy solid for ^Bragg = 0 , and dashed for 
negative values (gieen to blue color shades). 


and the conesponding duectional spread, 

ae= {2[l-(<7icose-hZ>isin0)/£-(/p)]}^ (m.74) 

O 0 ranges fiom 0 for miiduectional waves, to 81 degiees for isotiopic waves. Offshore 
propagating waves (it < 0 < 27t) are excluded m the analysis because the predicted back- 
scattering is weak, and reflection fiom the beach (Elgar, Herbers & Guza, 1994), not rep¬ 
resented in the model, is apparent in the 8 m data (figure 29b). 

The model without Bragg scattering predicts the expected timimg of 0 towards 
the shore-nonnal direction, caused by refiaction (figme 29a). The introduction of Bragg 
scattering shifts the mean wave direction by an additional 1 to 10 degrees to the north, 
because the bottom spectnim is not isotropic (figme 26). This effect is strongest for the 
hindcast which uses the bottom spectrmn Ff with a larger variance. This small shift is not 
evident in the observations, suggesting that other processes, not represented in the model, 
may be important or the orientation of the bathymetric featmes in figure (26a) may not 
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Figure 29. Diiectional properties across the shelf 

Measiued (•) and predicted (-, with Bragg scattering using -, with Bragg 

scattering using ( 1 ); ., without Bragg scatteiing) variations of 0 (a) and ae (b) 

at the peak frequency fp. Results are shown for October 20, 2000, at 2100 GMT, along a 
cross-shelf transect (c) extending from the 8 m depth anay to deep water offshore of NDBC 
buoy 44014 . A maximum value of k/I = 5 was used m the scattering calculations. 
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be representative of other parts of the shelf, but it is unlikely that the statistical nature of 
these features evolved significantly between 1994 when the wave data was collected and 
1999 when the bottom bathymetry was surveyed, even if individual features may have 
moved. Besides, this effect is still apparent for waves observed in 1999 (see § V.D.l). The 
detailed directional spectra, shown in figure 30 demonstrate that rather than shifting the 
entire spectrum, Bragg scattering skews the directional spectrum to the north (figure 30c,d) 
by preferentially scattering waves that propagate in directions parallel to the crests of the 
larger bedforms (i.e. waves from the north-east, figures 26, 26). 

Bragg scattering strongly affects the directional spread, causing a gradual increase 
of Oe across the shelf (figure 29a), that partly balances the reduction of the directional 
spread of the incident waves caused by refraction. Results based on bottom elevation spec¬ 
tra Ff and are qualitatively similar but the increase in directional spread is much larger 
for the more ‘energetic’ bottom spectrum F^ (about a factor 2.5) than predicted for F^ (a 
factor 1.6). On the inner shelf, in 8 m depth, the observed Oe value of 14° is a factor 2 
larger than the model prediction without Bragg scattering (7°, figure 29b), but falls in the 
range of model results with the source term Ssragg based on bottom spectra F^ (18°) and 
F^ ( 12 °). 

The cut-off value {k/l)^^ of the ratio between surface and bottom wavenumbers 
was varied from 0.5 (no scattering) to 5, in order to examine the importance of different 
bottom topography scales in the scattering process (figure 31). Increasing {k //) 

(no scattering) to 1 (maximum scattering angle 0max = 60°) does not change significantly 
the directional properties of the waves. These interactions, involving bottom components 
with wavelengths smaller than the surface wavelength, are weak because of the sharp roll¬ 
off of the bottom spectral levels at high wavenumbers. At the other end of the spectrum, 
results for {k/l)^^^ values of 4 and 5 are nearly identical, indicating that larger bottom 
features also do not significantly affect directional properties. Although the bottom spectral 
levels are relatively high at these small values of I, the angular separation of the interacting 
wave components is small (11.5° for k/l = 5) and thus the energy transfers do not strongly 
modify a directional spectrum that is already broad. A range of interactions involving 
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(d) model output, 8 m depth, with ^Biagg 
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Figiue 30. Realistic test of the Bragg scattering soiuce tenn 
Observed wave spectia in (77 ) 49 and (Zt) 8 m depth, and predicted wave spectra in 8 m 
depth, (c ) without Bragg scatteiing and (d ) with Bragg scattei'ing, based on the bottom 
spectnun Ff (figiue 2677 ). Note that waves co inin g from the west in panel b aie probably 
reflections fr om shore (at 1 km of the 8 m site). 
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Figme 31. Sensitivity to the wavenumber cut-off 
Predicted variations of {a) 0 {fp) and (Z») ae {fp) across the shelf for six values of the cut-off 
parameter -, 5; O, 4; x, 3; A, 2; -I-, 1;., 0.5 (i.e. no Bragg scattering). 
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intermediate scale bottom components (A:// = 1 to 4) appears to dominate the scattering 
process (i.e. note the gradual shift of 0 and increase of Oe as increases from 1 to 

4 in figure 31). Predictions of Oe are insensitive to the bathymetry smoothing scale used in 
ray computations, with a typical 2 degrees difference between model runs using the original 
150 m resolution grid, and the predictions presented here using a 2 km smoothed grid. 

D. SUMMARY 

The energy balance equation for random surface gravity waves, including Bragg 
scattering (the lowest order resonant interactions between waves and bottom undulations), 
was re-derived for non-stationary conditions and multiple-scale bottom topography, com¬ 
bining Hasselmann’s (1962) perturbation expansion of the wave energy, with a ray approx¬ 
imation for medium variations. The bottom topography is decomposed in a large scale 
topography, responsible for wave refraction and shoaling, and random undulations with 
smaller wavelengths (of the order of the surface wavelength), that cause Bragg scattering. 
The effects of the large-scale and small-scale bottom slopes, surface non-linearity, wave 
non-stationarity and non-uniformity are represented by five small parameters, |3, T], £, a, 
and y, respectively. Using a ~ (3 ~ y ~ £^, a closure of the fourth order energy yields 

a spectral energy balance equation in which refraction, shoaling, and Bragg scattering pro¬ 
cesses are all of the same order £^. 

The stochastic scattering theory was reconciled with a spectral application of deter¬ 
ministic theories for waves propagating over sinusoidal bars (Davies & Heathershaw, 1984; 
Mei, 1985; Dalrymple & Kirby, 1986). Agreement of the theories in the asymptotic limit of 
a large number of bars supports the present derivation of the Bragg scattering source term 
(III.48) which is a factor 8 and 4 smaller than expressions given by Hasselmann (1966) 
and Long (1973), respectively. Analysis of the detuning of wave-bottom interactions by 
the large scale bottom slope (3 shows that the present theory is valid only for small values 
of |3/ cos^G, where 0 is the wave incidence direction relative to the bedform-normal. In 
the present application, wave-bottom interactions corresponding to 0 larger than a cut-off 
value 0max are neglected. 
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The effect of bottom scattering on swell propagation was illustrated with a hindcast 
for the North Carolina continental shelf using the numerical wave model CREST with high 
resolution bathymetry and an efficient implicit scheme to evaluate the bottom scattering 
contribution to the energy balance. Back-scattering of waves towards the open ocean was 
found to be negligible in this region. However, forward scattering causes a diffusion of 
wave energy about the mean direction that results in a dramatic increase of the directional 
spread of the wave spectrum on the inner shelf. This weak back-scattering and strong 
forward scattering is caused by the sharp roll-off of the bottom elevation spectrum at high 
wavenumbers. The directional broadening of the swell spectrum predicted in shallow water 
is qualitatively consistent with field measurements. 
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IV. OBSERVATIONS OF WAVE-GENERATED 

RIPPLES 

This chapter was submitted for publication, with part of chapter I, in the Journal 
of Geophysical Research (Ardhuin, Drake & Berbers, submitted manuscript). It describes 
observation of bedforms in relation to wave measurements. 

A. EXPERIMENT 

Wave, bedform and sediment data were collected off the coast of Duck, North Car¬ 
olina, from September to December 1999, as part of the SHOaling Waves Experiment 
(SHOWEX) a multi-institution effort to improve the understanding of the transformation 
of waves in shallow water. The experiment site is a wide sandy continental shelf, exposed 
to North Atlantic swells. The bottom slopes gently from 20 m depth, 5 km from shore, to 40 
m, over a distance of 60 km (figure 32), and is characterized by the presence of sand ridges 
with spacings and heights ranging from about 1-10 km and 1-10 m, respectively (Green, 
1986; Wright, 1995; figure 25). The larger-scale ridges are oriented along a south-north 
axis (figure 32). 

Six Datawell Directional Waverider buoys, named XI to X 6 , were deployed from 
September 13 to December 13, 1999, along a cross-shelf transect. Buoys XI through X5 
span the shelf from 21 m to 39 m depth, and buoy X 6 was located on the shelf break in 193 
m depth. Wave frequency spectra and directional moments (mean propagation direction and 
directional spread) as functions of frequency were computed internally by the buoys at 30 
minute intervals, and transmitted continuously to shore by buoys XI to X4 using HE radio. 
Buoys X5 and X 6 , situated farther offshore, were equipped with internal data loggers. The 
significant wave height // 1/3 (four times the root-mean-square surface displacement) and 
spectral peak frequency fp observed at XI and X3 are shown in figure 33 for energetic 
periods preceding three cruises when side-scan sonar images of the bottom were acquired. 
Eigure 33 also displays the mean wave direction 0£ (/), defined as the direction of the first 
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Figure 32. Topography and SHOWEX sidescan sonar surveys 
Squares numbered XI to X6 indicate the positions of Directional Waverider Buoys 
deployed during SHOWEX, the large filled circle indicates the National Data Buoy 
Center buoy 44014, and the + symbol indicates the Field Research Facility 8 m depth 
pressure sensor array (8M). Solid lines represent ship tracks for the sidescan sonar 
snrveys, and triangles indicate the locations where surficial sediments were sampled. 


96 











a 



Figure 33. SHOWEX observed wave conditions and ripples 
(a) Evolution of significant wave height i/ 1 / 3 , (b) mean wave direction at the peak 
frequency 6e (fp), and (c) peak frequency fp. Observed are shown at buoys XI (solid) 
and X3 (dotted), spanning the inner shelf. Before survey Si (conducted during the 
buoys deployment), data from 8 M and 44014 are shown. During Hurricane Dennis 
&E (fp) is estimated to be 0-10° larger than 6e (fp) measured at 8 M, owing to refrac¬ 
tion. For each survey the average of crest-normal ripple directions 6r observed in the 
vicinity of XI is indicated with a solid horizontal bar. The dashed bar indicates more 
southerly average Or of ripples observed close to X3 during S 2 . 
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The large, long-period waves (peak frequencies ranging from 0.06 to 0.08 Hz) originated 
from three category 4 hurricanes: Dennis looped offshore of the buoy transect (August 29 - 
September 5) before coming ashore farther south, Floyd swept through the region, with the 
eye passing a few tens of kilometers inland on September 15, and Gert remained more than 
1200 km away from the experiment site, sending swell along the buoy transect (September 
20 to 23). The month of December was marked by a strong northeaster (December 1 and 2). 
The observed strong attenuation of wave energy between X3 and XI, particularly for low 
frequency waves when Hi/^, is in the range 1-2 m, is consistent with previous observations 
and the large bottom roughness predicted for wave-formed ripples (Herbers etal, 2000; 
Ardhuin et al, 2001). 

Bedforms were surveyed with an EG&G model DF1000 side-scan sonar, towed 
from the R.V. Cape Hatteras at night during cruises for buoy deployment, maintenance, 
and recovery. The first survey (SI) took place between the passages of hurricanes Dennis 
and Floyd from September 10 to 13. A second survey (S2) was conducted two weeks 
later (September 24-26), after swells from Gert and a coincident Northeasterly wind sea, 
had subsided. Repeating most of the tracks from survey SI, the objective of S2 was an 
assessment of changes in bedforms after hurricanes Floyd and Gert. The final survey (S3) 
took place as buoys were being recovered (December 12-15), after a series of northeasters. 
The same tracks were repeated again, and a few additional regions were surveyed (all ship 
survey tracks for S3 are shown in figure 32). 

The dual frequency sonar (100 and 500 kHz), was towed 2-3 m above the bottom 
in the shallower regions and 10-15 m in regions deeper than 30 m, at a speed of about 
2.3 m s^^ The 100 kHz frequency was used throughout for larger area coverage, at the 
expense of poorer resolution. All back-scatter sonograms, covering a 100-m wide swath 
along a total track length of 420 km were both recorded in a digital format and printed 
on paper scrolls. In selected regions where ripples were clearly visible, the digital data 
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Figiue 34. Example gjaiii size distiibutions 

A typical distribution (sobd thick), and the distributions with the finest (asterisks) and 
coarsest (triangles) sediments are shown. 

was transfonned into GEO-TIFF mosaics, fully collected for slant range and speed, with 
a resolution of 5 or 10 pixels per meter. Siuficial sediments were collected with a Shipek 
grab (46 samples), and by divers (5 samples), at locations indicated in figiue 32. 

B. SURFICIAL SEDIMENTS 

The continental shelf between Cape Hatteras and the entrance to the Chesapeake 
Bay (Cape Hemy) is believed to be covered by fin e quartz sand (Milliman et a I., 1972; 
Swift and Sears, 1974), with occasional coarser sands and gravel in fluvial paleochamiels 
flooded by Holocene sea-level rise (Swift et al., 1972). The Albemarle river paleochannel 
comcides with the buoy transect (figiue 32). Very fine sands with a silt and clay content of 
10 to 26% are foiuid between 8 and 18 m depth (Field et al., 1979; Wright, 1993; Madsen 
et oL, 1993), then origin can be traced to old lagoonal deposits (Wright, 1995). 

The samples gathered dining the three cruises confirm this pattern. Most samples 
consisted of fine sand with a median diameter D^q of about 0.2 rum (figiue 34). Samples 
occasionally contamed shell hash, particularly ar oimd XI (figiue 32). Grain sizes at neigh- 
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boring sample sites can be significantly different (figure 35a, suggesting strong variability 
of sediment properties on scales less than a few kilometers, in particular on the inner shelf 
between XI and X2. For example, the finest sand (D^q = 0.09 mm, asterisks in figure 
34) was found offshore of the Field Research Facility pier in 12 m depth, and the coarsest 
sample, a mixture of sand and small gravel with D^q = 4 mm (triangles in figure 34), was 
found in the vicinity of XI, in 19 m depth. All samples taken near X2 were fine-medium 
sand, except for one coarse sample collected by divers within 300 m of X2 that contained 
pebbles several centimeters in diameter. 

The variability of sediment properties also is evident in the side-scan sonograms. 
For the range of grain sizes found in the samples, a stronger reflectivity (dark shades in 
the sonograms) usually corresponds to coarser sands. A large portion of the regions sur¬ 
veyed between XI and X2 show alternating light and dark bands (the coarsest sand found 
in the samples came from one of these dark bands), a few hundred meters wide, with sharp 
transitions (e.g. figure 37a, discussed below). Offshore of X2 some dark bands are still 
present, although they contrast less with their lighter surroundings. This alternating pattern 
is similar to the one found by Green (1986) in a side-scan sonar survey conducted in 1984, 
6 km south of XL Green (1986) identified dark bands lying in the troughs between ridges 
oriented north-south (with approximate heights and spacings of 5 and 500 m, respectively) 
as the pre-Holocene ‘basement’ underlying the finer Holocene sand ridges. In the surveys 
presented here the dark regions also correspond to relatively low-lying areas, at least in¬ 
shore of X3 where high-resolution multibeam sonar surveys were conducted in December 
1999 (submitted manuscript). 

C. SIDESCAN SONAR IMAGES OE BEDEORMS 

In addition to qualitative information on sediment grain sizes, the sidescan sonar 
images document the presence and characteristics of bedform patterns. A qualitative in¬ 
spection of the entire sonar data set for the three surveys (figure 35c-e) reveals the predom¬ 
inance of long-crested ripples with wavelengths between 0.4 and 3 m, and crests approxi¬ 
mately parallel to the coastline, consistent with the scales and orientation of vortex ripples 
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Figme 35. Spatial variations of and ripple coverage 
(a) Geographical distribution of the median grain size (b) depth profile along the 
instnunented transect, and (c,d,e) fraction of sonar smveyed regions visibly covered by 
ripples for siuveys SI, S2 and S3, respectively. 
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formed by large low-frequency surface gravity waves propagating across the shelf. In cross¬ 
shore surveys ripple marks are generally faint or absent, probably because the alongshore 
orientation of the sonar beams (perpendicular to the ship track) was at a grazing angle rela¬ 
tive to the ripple crests, usually oriented alongshore. Images from part of the surveys were 
blurred by excessive wave-induced towfish motion. In regions where ripples were not de¬ 
tected, the bottom may be smooth, as a result of the shallow nature of the sediment cover or 
some other unknown phenomenon, or the bottom may be covered with ripples that are too 
small to be resolved. For example at a site in 12 m depth high-resolution mapping instru¬ 
ments located on a bottom-mounted frame show the presence of short wavelength ( 10-20 
cm) ripples (T. Stanton, Naval Postgraduate School, personal communication, 2000) in 
fine sediments (D 50 = 0.09 mm, asterisks in figure 34) that were not resolved in the sides- 
can sonar surveys. The absence of ripples in alongshore ship tracks usually correspond to 
regions with finer surficial sediments. Assuming that vortex ripples, for which the wave¬ 
length scales with the near-bed wave orbital diameter, are formed by wave action during 
a storm and become relic when the near-bed velocity decreases, ripples generated on finer 
sand are ‘frozen’ later, when the waves are smaller, and thus are expected to have shorter 
wavelengths (i.e. not resolved by the sonar) than ripples in coarse sand. A quantitative 
analysis was performed for surveyed regions with clearly visible bedforms and relatively 
uniform bedform geometry and sound reflectivity. For each region a representative 100 m 
by 100 m image was processed. In some cases with faint ripples a smaller image with clear 
ripple crests was analyzed. Although spectral analysis is well suited to characterize reg¬ 
ular ripples, bidimensional variance spectra computed from sonar images exhibit multiple 
peaks in many cases, either as the result of excessive sidescan towfish motions (along-track 
modulations of ripple patterns on the sonar images), or because of the presence of many 
defects in the ripple patterns. Therefore a different procedure was used that is less sensitive 
to ripple defects and image artifacts (figure 36). For each image, after subtracting a piece- 
wise bilinear fit to the pixel values, we computed the zero-crossing contours of the image 
intensity, and formed a histogram of the contour length as a function of its orientation, at 
one degree intervals, using the trigonometric (right-handed) convention. 
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Figure 36. Extraction of ripple parameters from sidescan images 
(a, b, c) example images from the sites where samples 111, 114 and 207 were collected 
(see table I), (d, e, f) corresponding zero-crossing contours with an indication of 
the mean crest orientation 6c (thick dashed line), and (g,h,i), histogram of contour 
direction with htted cosine distribution (thick dashed line). 
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A function p{Q) = ao + a 2 COS [2(0 — 0^)] was fitted to this histogram in a least- 
squares sense, and 0^ was interpreted as the mean orientation of the ripple crests. The ratio 
a 2 /aQ gives an indication of the presence, long-crestedness, and regularity of ripples. For 
high contrast ripples with long crests and few pattern defects ^ 2/^0 > 1 (e.g. figure 36a,b). 
For faint ripples or ripples with short or brick-patterned crests ^ 2/^0 < 0.5. An example 
of marginally detectable and non-uniform ripples, with = 0 . 8 , in shown in figure 

36c. Once 0c was determined, standard one-dimensional Fourier spectra were computed 
for east-west image lines (usually within 30® of the crest-normal direction). The spectra 
were then averaged over all lines, and the peak wavenumber kp was determined. The 
average wavelength of the ripples \ was then taken to be 27tsin0c/kp, and the crest-normal 
direction, using the (left-handed) nautical convention, was given by 0 ^ = 180(1 — Oc/tt). 
In cases with ^ 2/^0 < 0.4, this procedure failed to determine 0^ and \ unambiguously, and 
the bed was presumed featureless. 

In most images with visible ripples the dominant bedforms are regular long-crested 
ripples with an apparently sinusoidal profile (figure 36). Spatial variations observed in the 
backscatter sonar intensity and ripple patterns are the result of non-uniformities in sediment 
properties (e.g. variations in grain size (figure 37a) or small-scale (unresolved) ripples or 
both. A few regions reveal intricate bedform patterns with second harmonic ripples (figure 
37b) that may be the result of decreasing wave forcing conditions. Earlier observations 
and numerical modeling studies have shown that when the mean wave direction is con¬ 
stant, ripple wavelengths gradually increase with increasing wave orbital diameter J 1 / 3 , 
but the converse is not true: when decreases the ripple wavelengths typically remain 
unchanged. For decreasing J 1 / 3 , if the Shields number threshold for sediment motion is still 
exceeded, secondary ripple crests may appear in the troughs, reducing the wavelength by 
half (Forel, 1894; Traykovski et al., 1999; Andersen, 1999). Other ripple patterns include 
superpositions of two ripple systems (e.g. figure 37c,d) similar to previous observations 
(e.g. Forel, 1894; Swift etal., 1972). As discussed below, these patterns are likely caused 
by successive wave events with different wave directions. The wavelengths in some of 
these ripple systems can be as large as 4 m, of the order of the wavelengths of long-crested 
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Figure 37. Examples of unusual ripple patterns 
(a) sharp transition in sediment properties, (b) second harmonic ripples, and (c and 
d) superposition of ripple patterns (see text). 
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Figiue 38. Time evolution of ripple directions and wavelengths 
(a) for survey SI, (b) S2, and (c) S3. The thick bars indicate the crest-noiinal ripple di¬ 
rection 0;. and the r ipple wavelength X with a scale given in the upper right corner of each 
panel. Each bar is drawn in the middle of the survey segment for which the ripple param¬ 
eters are representative. The colors along the survey line indicate the average back-scatter 
strength of the sidescan sonar images: red is strong (generally coarse sand) and blue is 
weak (generally very fine sand). Siu'veyed regions withoirt a bar indicate the absence of 
r ipples or a faihue of the analysis procediue to determine 0^ and X. 

featmes observed in this region and called ‘rnegaripples’ by Green (1986). The variety of 
observed bedform patterns imderscores the important effects of past bedform evohrtion on 
then present state in the relatively mild forcing conditions prevalent dining these sirrveys. 

Temporal changes in ripple direction 0^ and wavelength X are examined in figine 
38 irsing repeated sinwey track lines close to XL All smveys reveal widespread ripple 
coverage with typical wavelengths X of 0.5-2 m, and ripple crest-normal directions 0^ 
between 60 and 120° (shore-normal is 70°). Changes observed between surveys SI and S2 
(figiue 38a,b) are small, with a slight shift of the average value of 0r from 91° to 85°. A 
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larger change in ripple orientation is observed in S3 with an average 0;- of 13° (figure 38c). 

This observed change in bedforms is related to a change in wave conditions, from 
hurricane-generated waves from east to south-east propagation directions in September, to 
northeasterly waves in December. The average values of 0^ in surveys SI, S2 and S3 (indi¬ 
cated by thick horizontal bars in figure 33) correspond to the wave directions of preceding 
large wave events on September 4-5 (Hurricane Dennis), September (22-23) (Hurricane 
Gert), and December 1-2 (northeaster), suggesting that the observed ripples were formed 
in these periods. The large oblique ripples observed in S2 (X = 2.0 m, 0^ = 105") and 
S3 (X = 2.0 m, 0^ = 132") close to the northern end of the surveyed region are notable 
exceptions. The corresponding sidescan image of the anomalous ripples in S2, shown in 
figure 37c, indicates that the estimated X and 0^ are averages of two superposed ripple sys¬ 
tems: X = 2.0 m, 0^ = 120", probably generated during the passage of Hurricane Floyd, and 
X = 0.8 m, 0^ = 85" probably generated during the arrival of swell from Hurricane Gert. 
Within 200 m of these ripples, a similar but smaller patch of large ripples was observed in 
survey S3 (figure 38c) with a more oblique angle 0^ = 132" that does not match the wave 
directions of any energetic wave event observed during the three months between surveys 
S2 and S3 that would be capable of moving bottom sediments in that region, suggesting 
that these ripples also were formed during Floyd and persisted throughout the experiment. 

Spatial variations in 0^ and X across the shelf during survey S2 are shown in fig¬ 
ure 39. Observations on the middle and outer shelf show more southerly 0^ angles than 
observations on the inner shelf (figures 39b,c and 39a). Some spatial variation of ripple 
directions is expected associated with the refraction of the swell that probably generated 
the bedforms. However, observed differences in wave direction fig between X3 and XI are 
small compared with the 15" difference between the average 0^ around X3 and XL Tempo¬ 
ral changes of the wave forcing conditions also may have contributed to the 0^ variations. 
In deeper regions where the wave motion is more attenuated over the water column, ripples 
may have become relic at an earlier stage when the swell was more energetic and arrived 
from a more southerly angle (figure 33a,b). Indeed the average 0^ values of 100" and 85" 
near X3 and XI, respectively (dashed and solid horizontal bars in figure 33b), are close to 
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Figiue 39. Cross-shelf dishibution of ripple dir ections and wavelengths 
Same as figme 38, for siuvey S2 only, at different locations across the shelf. 

mean wave duections on September 22 and 23, respectively. Large variations of 0^ and X 
also are observed over much shorter distances, with estimates varying by as mitch as 30° 
and a factor 2, respectively, within 1 km. These differences are cairsed primarily by the 
superposition of ripple patterns (figures 39b and 39c), some of which are barely resolved 
in oin images. The relation between observed ripple geometry and changing wave forcing 
conditions is analyzed in the next section. 

D. EVALUATION OF RIPPLE PARAMETERIZATIONS 

Bedfornis observed in the sonar images irsirally are not related to the instantaneous 
wave forcing (sinveys were performed hr relatively calm conditions), but are the result of 
earher wave events that were sirfficiently energetic to move the sediments. Here, the ripple 
properties observed in sinveys S2 and S3 are compared with the wave forcing history es¬ 
timated at the nearest buoy. Neglecting variations ui the wave forcuig conditions over the 
1-5 km distance separating the locations of the wave and r ipple measiuements, we com¬ 
puted a ‘significant’ Shields number nri /3 from (1.4), using a range of D^q values for the 
grain size, and ;/i /3 defined as two times the root-mean-square velocity near the sea bed, 
determined fiom the sinface elevation frequency spectrum with linear wave theory. The 
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skin friction factor was obtained with Grant and Madsen’s (1979, 1982) parameteriza¬ 
tion, based on a linear profile of the eddy viscosity (see also the review by Tolman, 1994). 
We assumed a threshold of ripple formation \|/i/3 = \|/c (Traykovski et al, 1999), where t|/c 
is the threshold of sediment motion under sinusoidal waves. Madsen and Grant (1976) give 
a review of values of determined experimentally for well-sorted sand. Here we use a 
piece-wise fit to the experimental data of Wallbridge etal. (1999) for mixed sands under 
combined wave-current flow, that we apply to the median grain size, 

vitc = 0.1 exp [(S*-2) In (0.35)/10] for 2 < S* < 12, (IV.2) 

S* = D5o[(^-l)gZ)5o]^/(4v), (IV.3) 

where V is the kinematic viscosity of water. These values of t|/c are slightly larger than 
those given by Shields (1936) or Soulsby and Whitehouse (1997) for unidirectional and 
combined wave-current flows (about 30 % larger in the range 3 < 5* < 12), but follow the 
same trend, decreasing with increasing grain size. Based on these assumptions, all but the 
coarsest surficial sediments have moved on 22 September, and during the 1-2 December 
northeaster, after which ripples became relic and remained unchanged until surveys S2 and 
S3, respectively (figure 40a). 

The interpretation of survey S2 is complicated by cross seas observed on September 
22 (when the wave height was maximum). Different definitions of a representative wave 
forcing direction yield different ripple predictions depending on the relative importance 
they give to motions with higher and lower frequencies. Here the bulk mean directions of 
the near-bed orbital velocity Bu and displacement 0^ are defined as the directions of the first 
spectral moment vectors are defined as 

{au,bu) = [ . [ (cos0,sin0) £(/,0)d0d/, (IV.4) 

J smh^ [kH) Jo 

f 1 /■2;r 

{ad,bd) = / ., 2 nzj\ (cos0,sm0) £(/,0)ded/, (IV.5) 

J smh^ [kH) Jo 

where the angle 0 is defined with the nautical convention. The superposition of a local 
wind sea traveling at right angles to swell from Hurricane Gert caused to veer to the 
north from 85 to 50® and then back to the east, whereas 0^ shifted gradually from 95 to 70® 
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Figiue 40. Wave forcing conditions dining ripple foiniation 
(a) y 1/3 /v|/c computed at X1 using the flill wave spectnim and various values of the median 
gr ain size D^o- (b) Coiiesponding values of the mean dir ection of the orbital displacements 
at the top of the wave bottom boimdary layer Qj. For the complex September 22-23 wave 
conditions the mean dir ection of the near-bed orbital velocity 0„ is also indicated, as well 
as Oj estimates for swell (0.05-0.11 Hz) and wind sea (0.11-0.25 Hz) only, (c) Cone- 
sponding values of predicted ripple wavelength 0.1 for vorfex ripples, where dij^ is 
the significant diameter of the orbital displacement at the top of the wave bottom boimdary 
layer. The gray band in (b,c) represent the mean value of the crest-normal r ipple direction 
0,- and wavelength X, respectively, plus or minus one standard deviation, observed in the 
vicinity of XI (figme 38b,c). 
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(figure 40b). Although the high frequeney seas eontributed signifieantly to the near-bed 
veloeity varianee and Shields number \|/i /3 (equal eontributions from the 0.05-0.11 and 
0.11-0.25 Hz frequeney ranges when is maximum), the orbital displaeement diree- 
tion 0^ remained aligned with the lower frequeney swell. Observed ripple erest-normal 
direetions 0^ are eloser to 0j than 0„ at possible freezing times (\[ti /3 ~ \[tc, figure 40a). 
These observations suggest that the mean direetion of the orbital displaeement 0^, whieh 
is less influeneed by high frequeney (wind sea) motion, may deseribe better the relation 
between wave direetion and ripple orientation (figure 40b). The V|/i /3 = \\fc ripple freezing 
threshold yields values of 0^ that are still 5-10® to the south of 0^. Several faetors may 
eontribute to this differenee: a reorientation of the ripple pattern may require a larger fore- 
ing (tl/ 1 / 3 ) than the initiation of sediment motion (t|/c) (Werner & Koeurek, 1997), the time 
of adjustment of the bedforms in eonditions elose to the threshold may be larger than the 
time seale of the wave height deerease, 0 ^ may not be simply related to Qj, or the values 
of \|/c are not well predieted by (IV.2)-(IV.3). Sueh ‘early freezing’ of ripples in eonditions 
of deereasing foreing also was observed by Traykovski eta/.(1999) with frozen bedform 
patterns oeeurring for \|/i /3 as large as 2\|/c . Although the exaet threshold for the eessa- 
tion of ripples evolution is uneertain, the wavelength X are generally elose to when 

¥c < Vi /3 < 3\|/c for representative grain sizes (figure 40a,e), suggesting that the observed 
ripples are probably of the ‘vortex ripple’ type, similar to those observed by Traykovski 
et a/.(1999). 

Analysis of bedforms at sites that eoineide with the loeation of sediment samples 
shows that the ratio X/D^q varies from 300 to 8100 (table I), extending the range of previ¬ 
ously observed values (Wiberg & Harris, 1994), in partieular for vortex ripples. Crude 
estimates of aetive ripple foreing eonditions ean be obtained from the preeeding wave 
reeords by determining the most reeent time when ripples and wave erests were aligned 
(|0rf-e,| < 5®), while surfieial sediments were still in motion (assuming t|/i /3 > 0 . 8 t|/c). 
This analysis gives values of between 0.9 and 3.5, and values of X/dij^ from 0.5 to 

0.76 (table I). Although these estimates have a large uneertainty and should be eonsidered 
with eaution, they are eonsistent with aeeurate measurements by Traykovski eta/.(1999) 
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Sample number 

111 

no 

114 

11 

207 

212 

206 

205 

Nearest buoy 

XI 

XI 

XI 

XI 

X3 

X3 

X3 

X3 

Depth (m) 

22.3 

19.7 

20.7 

21.2 

22.7 

27.6 

26.3 

27.4 

D 50 (mm) 

4.71 

0.98 

0.32 

0.15 

0.28 

0.19 

0.17 

0.15 

D 85 (mm) 

7.30 

2.29 

0.48 

0.28 

0.68 

0.28 

0.25 

0.24 

Br 

94° 

85° 

00 

0 

0 

00 

00 

0 

101° 

99° 

n.a. 

n.a. 

X (m) 

1.25 

1.20 

0.77 

1.22 

1.30 

1.37 

n.a. 

n.a. 

di/3 (m) 

2.55 

1.85 

1.34 

1.60 

2.73 

2.34 

n.a. 

n.a. 

mi/ 3 (ms~i) 

0.58 

0.48 

0.37 

0.42 

0.67 

0.55 

n.a. 

n.a. 

¥1/3 M 

0.90 

2.00 

1.33 

3.54 

3.4 

3.0 

n.a. 

n.a. 

Ws/ui/3 

0.50 

0.25 

0.10 

0.03 

0.08 

0.04 

n.a. 

n.a. 

V^l/3 

0.49 

0.65 

0.58 

0.76 

0.48 

0.59 

n.a. 

n.a. 

^1/3/-D50 

540 

1900 

4200 

10700 

9800 

12300 

n.a. 

n.a. 

V-Dso 

270 

1200 

2400 

8133 

4600 

7200 

n.a. 

n.a. 


Table 1. Ripples at sites of sediment sample eolleetion 
The foreing eonditions (orbital diameter and eorresponding veloeity and Shields num¬ 
ber) at the time of ripple ‘freezing’ were erudely determined by matehing the bottom orbital 
displaeement direetion 0 ^ at the nearest buoy to the direetion normal to the ripple erests 0 ^ 
(see text). 


who observed that X is elose to O. 76 J 1/3 under moderate foreing eonditions. These esti¬ 
mates are also qualitatively eonsistent with the values ofX/d predicted in numerical simu¬ 
lations with sinusoidal waves (Andersen, 1999). However there is no clear increase ofX/d 
with the normalized fall velocity Ws/u. 

Estimates of ripple direction 0^ and wavelength X in surveys S2 and S3 are sum¬ 
marized in figure 41 for the region around XI. The observations in S2 show a clear trend 
of increasing wavelength corresponding to more southerly angles (figure 41a). This trend 
is consistent with the freezing of the largest ripples shortly after the peak of the event 
(September 22), and a later freezing in finer sediments where the ripple patterns gradually 
adjusted to the more northerly angles, and small orbital diameters of the waves observed 
on September 23. Observations in S3 (figure 41b) do not indicate a clear relation between 
Br and X, possibly because the wave direction remained steady when the storm subsided 
on December 2. However, with the exception of apparent ‘Floyd survivors’ (0;- = 132°, 
X = 2 m), the longer ripples {X> 1.5 m) in S3 consistently have crest-normal directions at 
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Figiue 41. Ripple wavelength X versus crest-noiinal ripple direction 0^ 

X versus 0,- for ripples observed on the inner shelf in (a) siuvey S2 (figiue 38b), and (b) 
sirrvey S3 (figiue 38c). 
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large northerly angles (60® < 0^ < 75®), suggesting they became relic on December 1 and 
early December 2 when wave directions were similar, and the orbital diameter dij^, was 
maximum. 

E. SUMMARY AND DISCUSSION 

Widespread formation and evolution of sand ripples on the North Carolina conti¬ 
nental shelf, observed in side-scan sonar surveys, were examined in relation to time-series 
of directional wave measurements. Sediment samples confirm the ubiquitous presence of 
very fine to fine quartz-dominated sand, with significant variability in the grain sizes over 
distances less than 1 km, in particular in regions shallower than 25 m. Side-scan sonar sur¬ 
veys reveal the presence of ripples in most surveyed areas, with wavelengths between 0.5 
and 3 m, and crest-normal directions within 30 degrees of the west-east axis. The apparent 
absence of any bedforms in some regions with finer sediments may be caused by their close 
spacing (less than 40 cm) not resolved by the sonar images. 

The ripple crest-normal directions 0^ are aligned approximately with the direction 
of the near-bed orbital displacement 0^ when the ripples became relic. Although the precise 
time of cessation of sediment motion is uncertain, estimated ripple wavelengths X and 
forcing conditions (the orbital diameter J 1 / 3 ) are consistent with previous observations of 
vortex ripples ()i/Ji /3 ~ 0.5-0.76), even for large values of J 1 / 3 /D 50 and X/D^q. 

The present observations extend the range of previous field observations of vortex 
ripples (e.g. Traykovski etal, 1999) to larger adimensional orbital diameters J 1 / 3 /D 50 , 
but equally moderate Shields numbers 1 < \[ti/ 3 /\|/c < 3, for which similar vortex ripples 
(X/d ^ 0.8 with X/D^q up to 2 X 10^) have been observed in the laboratory (Southard et al, 
1990). Both the present observations and these earlier laboratory and field studies do not 
support parameterizations of the ripples based on the grain size only (e.g. Wiberg & Har¬ 
ris, 1994), but rather support a need for more complete parameterization of the physical 
processes responsible for ripple formation, including a dependence of ripples character¬ 
istics on the Shields number and possibly the sediment fall velocity (e.g. Nielsen, 1981; 
Andersen & Fredspe, 1999). 
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The widespread presence of vortex ripples on the continental shelf lends support 
for a dominant role of these rough bedforms in the strong attenuation of swells observed 
in this and other coastal experiments, during periods of moderate to strong wave forcing 
conditions (Hasselmann et al, 1973; Young and Gorman, 1995; Herbers et al, 2000), when 
ripples should be generated (Ardhuin et al, 2001). 
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V. VALIDATION OF THE ENERGY BALANCE 

EQUATION 


A. INTRODUCTION 

As described in chapter I many processes affect the evolution of gravity waves 
propagating at the surface of oceans, marginal seas or lakes, particularly in shallow water. 
Winds directly force the high frequency (‘wind sea’) part of the wave spectrum, actively 
generating waves with phase speeds slower than the wind speed, and a large fraction of this 
energy input is immediately lost through wave breaking (whitecaps). The low-frequency 
(‘swell’) part, with larger phase speeds, is influenced by local winds indirectly through non¬ 
linear wave-wave interactions. In water shallower than the surface wavelength additional 
wave-bottom interactions become important. As local winds dominate the wave climate in 
marginal seas such as the North Sea, investigations of the spectral energy balance of waves 
often have to deal with the full complexity of air-sea-bottom interactions (e.g. Bouws and 
Komen, 1983; Weber, 1988; Johnson and Kofoed-Hansen, 2000). However, outside the 
surf zone and in the absence of a wind-sea or strong currents, the steepness of swell is gen¬ 
erally too small to induce wave breaking, and the propagation distances across continental 
shelves are typically too short for significant nonlinear transfers of energy in the spectrum. 
Hence swell transformation across the shelf is a much simpler problem controlled primarily 
by wave-bottom interactions, and thus more tractable. Nevertheless, few studies of swell 
evolution in shallow water have been presented (e.g. Hasselmann et al, 1973; Young and 
Gorman, 1995; Herbers et al, 2000; chapter II). 

Swell is modified by a wide range of bottom topography and roughness scales, from 
the large scales (1-10 km) causing refraction and shoaling, and intermediate scales respon¬ 
sible for wave-bottom Bragg scattering, to small scales (0.1-10 m) enhancing the bottom 
roughness and wave energy dissipation, (figure 2). In earlier chapters a process-based ap¬ 
proach was used to estimate the dissipation of wave energy across the shelf (chapter II), 
related to the widely observed formation of sand ripples under waves (chapter IV), and to 
evaluate numerically the directional diffusion of waves, predicted by Bragg scattering the- 
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ory for random waves and irregular topography (ehapter III). Here we propose to eombine 
these proeesses in general swell eases, and test the elosure of the speetral energy balanee 
of swell with hindeasts of all swell-dominated periods observed during the experiments 
DUCK94 and SHOWEX. A eomprehensive statistieal analysis of the results is presented to 
evaluate different souree term formulations. This work will hopefully lead to more aceu- 
rate operational swell foreeasts on open oeean eoastlines, and is a first step towards elosing 
the energy balanee in more complex cases with wind forcing. An expanded version of the 
present chapter will be submitted as a two-part paper to the Journal of Physical Oceanog¬ 
raphy (co-authors: T. H. C. Herbers, P. F. lessen, and W. C. O’Reilly). 

We propose to test the following energy balance equation for the evolution of swell 

spectra 

—= 5'fric (k) + Ssragg (k) , (V. 1) 

where E is the wave energy spectral density, Sfhc and Ssragg are source terms representing 
bottom friction and wave-bottom Bragg scattering (III.48). Equation (V.l) incorporates the 
combined effects of all the scales of the bottom topography on a random wave field. The 
Eagrangian time derivative on the left hand side follows a wave component along its ray 
trajectory predicted by linear refraction theory (e.g. Mei, 1989). For Sfnc = 0 energy is con¬ 
served and (V.l) is identical to (III.47) derived in chapter III from the equations of motion 
assuming irrotational flow. Sfric , based on Tolman’s (1994) movable bed parameterization 
of wave-bedform coupling over a sandy bottom, is introduced heuristically in the energy 
balance as a sink term. A formal derivation of (V.l) with bottom friction may be performed 
by adding a rotational motion, following the approach of Weber (1991), to the irrotational 
flow determined in § III.A, but other source terms may arise as the result of the coupling 
of the rotational and irrotational flows. The bottom friction source term used here is further 
simplified by assuming a narrow wave spectrum (see for example Weber, 1991; Madsen, 
1994). 

We integrate (V.l) in time from deep water to the shore with the Coupled Rays with 
Eulerian Source Terms (CREST) wave model (chapter II), and the results are compared to 
wave measurements acquired during DUCK94 and SHOWEX on the North Carolina conti- 
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nental shelf. The wave data collected during the experiments, and the reduction of the data 
set to swell cases are described in § V.B, while the model set-up is detailed in § V.C. The 
results of hindcasts are given in § V.D, in the form of a statistical analysis of all observed 
swell events, followed by a discussion and summary in § V.E, where practical advice is 
given for swell forecasting in shallow areas. 

B. SWELL OBSERVATIONS DURING DUCK94 AND SHOWEX 
I. Instruments and data analysis 

During the DUCK94 (August-December 1994) and SHOWEX (September-December 
1999) experiments extensive wave measurements were collected across the wide shelf of 
the Mid-Atlantic Bight, in the region between Cape Hatteras and the entrance to the Chesa¬ 
peake Bay (figure 42). Bottom mounted pressure sensors were used during DUCK94 (Her- 
bers et ah, 2000; chapter II), and replaced by surface-following Datawell Directional Wa- 
verider buoys in SHOWEX. Eor the periods covered by both experiments wave data is 
also available from National Data Buoy Center (NDBC) 3-m pitch and roll discus buoy 
number 44014, in 49 m depth, and a coherent array of bottom pressure sensors and a Wa- 
verider buoy, in 8 and 15 m depth respectively, both maintained by the US Army Corps 
of Engineers Eield Research Eacility (ERE), in Duck, North Carolina. In addition to these 
instruments, the Diamond Shoals (DSEN7) and Chesapeake Eighthouse (CHEV2) C-MAN 
stations, operated by NDBC, were equipped with infrared laser wave gauges at the time 
of SHOWEX. The locations of the instruments are indicated in figure 42, and their basic 
characteristics are summarized in table II. 

The dataset includes both directional (buoys and coherent pressure array) and non- 
directional (single pressure, buoy, and laser gauges) measurements. Although some instru¬ 
ments are ideally equivalent (e.g. directional Waveriders and 3 m discus buoys), different 
sampling frequencies, record lengths, and response characteristics, introduce some vari¬ 
ations in quality. The model-data comparisons presented in § V.D are restricted to bulk 
spectral moments (significant wave height, mean direction and directional spread at 
the spectral peak, 0p and Oe p) that can be reliably estimated from the short records of 
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Figme 42. Bathymetiy and instnunent locations dining DUCK94 and SHOWEX. 
(see text and table II for instrament description) 
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Name 

water depth 

directional 

data availability 

Operated by 

8 M 

8.0 m 

yes 

02/1987 - present 

ERE 

WR(ERE) 

17.0 m 

no 

10/80-11/96 

ERE 

WR(ERE) 

17.0 m 

yes 

11/96 - present 

ERE 

44014 

49 m 

yes 

10/90 - present 

NDBC 

CHEV2 

15 m 

no 

9/84 - present 

NDBC 

DSEN7 

18 m 

no 

11/88 - present 

NDBC 

A 

12 m 

no 

DUCK94 until 17/11/94 

NFS 

B 

21 m 

no 

DUCK94 until 17/11/94 

NFS 

C 

26 m 

no 

DUCK94 

NFS 

D 

34 m 

no 

DUCK94 

NFS 

E 

35 m 

no 

DUCK94 

NFS 

E 

33 m 

no 

DUCK94 

NFS 

G 

46 m 

no 

DUCK94 

NFS 

H 

49 m 

no 

DUCK94 

NFS 

I 

87 m 

no 

DUCK94 

NFS 

XI 

21 m 

yes 

SHOWEX 

NFS 

X2 

24 m 

yes 

SHOWEX 

NFS 

X3 

26 m 

yes 

SHOWEX 

NFS 

X4 

33 m 

yes 

SHOWEX 

NFS 

X5 

39 m 

yes 

SHOWEX 

NFS 

X6 

193 m 

yes 

SHOWEX 

NFS 


Table 11. Wave-measuring instruments during DUCK94 and SHOWEX 
DUCK94 eovers 1/8/1994 - 30/11/1994, and SHOWEX spans 13/9/1999 - 13/12/1999. 
Oceasional instrument failure or maintenanee eausing loss of data is not indieated 


routine wave measurement systems. These parameters are determined from the wave fre- 
queney speetra E (/), and the standard lowest order Eourier moments of the direetional 
speetrum a\{f), b\{f), a 2 (/), and b 2 (/), computed for all directional instruments. Eor 
the high-resolution 8M array the full frequency-directional spectra E (/, 0) were estimated 
and offshore propagating waves were excluded from the directional moments. E (/, 0) 
spectra also were estimated from the directional moments at X6 and/or 44014 using the 
Maximum Entropy Method (Eygre & Krogstad, 1986), in order to initialize the model at 
the offshore boundary. All spectra and directional parameters were interpolated on a com¬ 
mon frequency grid, also shared by the numerical model CREST in the hindcasts described 
below. These parameters were averaged over 2-3 hour periods for the DUCK94 experi- 
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ment (variable length data records were collected every 3 hours), and 1 hour periods for 
SHOWEX (continuous data collection). The statistical parameters computed below from 
observations and model results are weighted accordingly, on an record length basis. 

2. Wave conditions 

The dominant wave events in the mid-Atlantic bight are generally the result of trop¬ 
ical storms and hurricanes in the summer-early fall, that follow a curved path to the north¬ 
west along the coast, or Northeaster storms that develop over North America and move 
offshore into the North Atlantic. 1994 was a moderate hurricane season, with only one ma¬ 
jor hurricane (Gordon, 15-19 November, figure 43) that remained for a long period in the 
Atlantic, south of Cape Hatteras, causing overwash of the barrier islands, and a maximum 
observed wave height of 9 m at 44014. Other significant events observed during DUCK94 
were Northeaster storms (figure 44). In contrast 1999 was a very active Hurricane season in 
the Atlantic. SHOWEX was delayed by a few days due to category 4 Hurricane Dennis (30 
August to 5 September) causing overwash in Hatteras Island. Instruments had just been 
deployed when category 4 Hurricane Eloyd made landfall south of Cape Hatteras, with 
maximum significant wave heights Hs of 6.75 m at 44014, 9 m at X 6 , and 12.5 m at NDBC 
buoy 41004 (34°30' N, located off Charleston, S.C., and not shown), and peak frequency 
/p = 0.11 Hz (at X 6 ). A maximum sustained wind speed t/ 19.5 = 34 m s~^ was recorded 
at the end of the ERE pier. Eloyd was immediately followed by Hurricane Gert (category 
4) that remained far offshore, sending large amplitude swell over the shelf (/p about 0.07 
Hz, Hs up to 3 m at X 6 , 3.4 m at 44014). and Hurricane Irene, that crossed the Elorida 
peninsula from the Gulf of Mexico into the Atlantic and passed 100 km offshore of Cape 
Hatteras, moving to the northeast. Hurricane Jose followed a track close to that of Gert 
(figure 43), sending only small amplitude swell to the instrument sites. SHOWEX was also 
marked by Northeaster storms, with a particularly severe storm on 1 December (figure 45). 
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Figiu e 43. Tracks of North Atlantic himicanes 
The dates indicate the position of the eye of the stonns (Gordon dining DUCK94, Floyd, 
Gei1, Irene, and Jose dining SFIOWEX), at 12:00 EST every day, after reaching the tropical 
depression stage. 
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Figme 44. Wind and wave conditions dming DUCK94. 

(a) Mean direction at the peak frequency fp, and (b) siguifrcant wave height, for buoy 44014 
and the 8M anay. Dark bands m (b) mdicate swell-douimated conditions, defrned here as 
periods when (c) the wind speed at 19.5 m is less than 60% of the wave phase speed at 
the peak frequency on the i nn er shelf, (d) Frequency spectmm (high densities m red, low 
densities in piuple, with a normalized logarithmic scale), and fp at buoy 44014. 
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Figiue 45. Wind and wave conditions dining SHOWEX. 
Same foimat as figiue 44. 
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3. Dataset reduction 

In order to restrict the present analysis to swell-dominated conditions in which the 
effects of wind and wave breaking can be neglected, the following criteria where applied 
to select hindcast periods: 

• a peak frequency fp less than 0.12 Hz (reduced to 0.10 Hz for DUCK94, to avoid 
potentially large errors in the depth correction of high frequency bottom pressure) 

• a maximum wind speed less than 0.6 times the linear wave phase speed at the peak 
frequency C(/p). 

The latter condition was chosen to exclude low-frequency wind-waves generated on the 
shelf in high wind conditions, such as the Hurricane Floyd landfall. The estimate of C {fp) 
are based on data from the Waverider buoy WR(FRF), on the inner shelf, and the wind 
speed is taken as the maximum of 1-hour averaged values [/ 19.5 measured at 19.5 m above 
sea level at the end of the FRF pier (close to the 8 M array) and t /5 measured at 5 m above 
sea level on board NDBC buoy number 44014. These instruments functioned with very few 
interruptions during both experiments. Wind speeds measured at the FRF pier and 44014 
are generally close, in spite of the difference in anemometer height, except for a time lag 
due to the motion of the weather systems and they are expected to bracket wind speeds at 5 
m height over the entire shelf. The wave phase speed estimated in 15 m depth give a lower 
bound for wave speeds on the shelf. Hence these criteria are expected to yield a conserva¬ 
tive selection of swell dominated conditions. The swell criteria were applied not only for 
the selected data record, but also during the preceding 3 hours, a period that corresponds 
to the propagation time of 0.08 Hz waves across the shelf. These swell-dominated time 
periods (indicated in black on figures 44b,45b) are generally associated with distant storms 
or the early arrival of low-frequency waves from approaching storms (figures 44d,45d). 

C. HINDCASTS ORGANIZATION 
I. The CREST wave model 

CREST is a hybrid Eulerian-Lagrangian spectral wave model described in chapter 
II, and § III.C, that uses finite frequency-direction bands, and an unstructured geographical 
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grid. The spectral energy balance equation is solved in its Lagrangian form (V.l), a simple 
one-dimensional advection equation, by advecting the k-space spectral densities of the 
wave energy along precomputed ray trajectories, from the model domain boundary to each 
point of the geographical grid. The source terms in (V.l) are determined in Eulerian form 
from the spectrum at each grid point and interpolated in space and directions on the rays to 
give Lagrangian source terms. For each grid point and finite band of frequency and arrival 
direction at that point (V.l) is integrated in time using a first order Euler scheme (§ II.A) 
with a fixed lO-minute timestep. 

2. Bathymetry and model grid 

In order to encompass the instrumented sites in both experiments, and allow for 
oblique swell arrivals, the model domain was set to cover the North Carolina and Virginia 
shelf region, extending 400 km between 34°30^ N, south of Cape Hatteras, and 38° N, at the 
Virginia-Maryland border on Assateague Island (figure 46). The National Ocean Service 
(NOS) database of depth soundings in that region generally has a resolution better than 200 
m, with finer details resolved around shoals dangerous for navigation. A large gap in this 
dataset, off Duck between X2 and X5, was filled with soundings acquired during instrument 
deployment, maintenance, and recovery, and dedicated cruises on board R/V Cape Hatteras, 
together with multibeam sonar surveys acquired on board the Canadian hydrographic vessel 
Frederick G. Creed (figures 11, 25). This composite dataset was gridded with 6” resolution 
in latitude and longitude (180 and 150 m, respectively). Before computing wave rays the 
bathymetry was smoothed with an isotropic tapered filter with a radius equal to 5 times 
the local wavelength (adapted to the wave frequency and local water depth), in order to 
remove smaller scale topographic features that are represented in a stochastic form in the 
wave-bottom Bragg scattering source term (see §III.B.4). The unstructured model grid 
was constructed starting from the instrument positions. Points were added along the 8, 15, 
20, 25, 30, 40, 60 ,100, and 1500 m isobaths, increasing spacings away from the coast 
and away from the instrumented transects. This ensemble of points was complemented by 
gradually introducing points at the centers of the corresponding Delaunay triangles, until 
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Figiue 46. Model giid. 

The grid points from which rays are computed and where the somce temis are evaluated 
are the nodes of the triangular' mesh. A linear interpolation is applied in each triangle to 
approximate the somce terms along the rays. The entire model domain is divided into 
sirbdomains, nimibered 1 throrrgh 9, separated by thicker lines. The locations of some 
insh'rrments are added for geographical reference, and a dotted line marks the 100 m depth 
contom. 
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all sides of these triangles were shorter than 2 to 20 km (depending on the water depth 
and the alongshore location). The model uses 29 frequency bands from 0.05 to 0.15 Hz, 
and 76 direction bands regularly spaced at 5-degree intervals. For each of these finite 
bands, bundles of rays were traced from all grid points. The grid was subdivided into nine 
subdomains (numbered 1 to 9 on figure 6) coupled at their mutual boundaries, where the 
ray computations are stopped, as described in § II.A.l. 

3. Boundary conditions 

At the external boundaries bordering domains 1 (land), 2 (offshore), 3 and 4 (north 
and south model limits) the following conditions were applied. The offshore frequency- 
direction wave spectrum is interpolated from the MEM-estimated spectra at X6 and 44014 
(back-refracted to deep water, assuming parallel bottom contours). To examine the sen¬ 
sitivity to offshore boundary conditions, results are also presented using only X6 data. 
Although the offshore conditions generally varied slowly on time scales of several hours, 
the deep water spectra are computed at 10-minute intervals, in order to match the model 
time step At, using a linear interpolation. Time lags based on the deep water group speed 
of linear waves are applied between the boundary grid points and the measurement loca¬ 
tion for each frequency-directional band, assuming uniform conditions along wave crests. 
This procedure was also used for the lateral boundary conditions (between subdomains 4 
and 8, 9 and 3, see figure 46), crudely accounting for refraction of the offshore spectrum 
by assuming parallel depth contours and applying Snel’s law. This treatment of the lat¬ 
eral boundaries does not represent accurately the propagation of waves across shelf regions 
outside the model domain, but avoids artificial shadow regions created by closed lateral 
boundaries (a problem encountered in chapter II). At the boundaries with domain 1 (land) 
a zero incoming flux was prescribed, corresponding to the absence of wave reflection from 
the beach and surf zone. 

4. Bottom sediments and small scale topography. 

Qualitative inspection of 20 core samples collected in 1997 on the inner shelf (Re¬ 
becca Beavers, Duke University, personal communication, 1999), and quantitative grain 
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size analysis of 50 samples gathered during SHOWEX (see also § IV.B), give a good de¬ 
scription of the surficial sediments in a narrow (20 km wide) region around the 8M-X5 
transect (figure 47). The median grain diameter D^q varies between 0.09 and 4 mm. The 
average value of (binned as a function of distance from shore) is 0.22 mm, slightly 
coarser than the average of D^q for the inner shelf (0.15 mm) used in chapter II. In order to 
extrapolate these observations away from the 8m-X5 transect a 5* order polynomial was 
fitted to the distribution of D^q as a function of the logarithm of the distance from shore 
(figure 47). This fitted distribution is used in model hindcasts, although very similar results 
obtained with a uniform value D^q = 0.22 mm suggest that the model is insensitive to the 
weak spatial variations in median grain size. 

The small-scale bottom topography is generally well resolved in the NOS bathy¬ 
metric database, for water depths between 10 and 30 m, and wavelengths larger than 300 
m. Estimates of bottom elevation spectra from a region with shoals on the inner shelf 
south of the instrument transects, and a deeper (20-40 m) mid-shelf region to the south¬ 
east of X2, give direction-integrated spectra with a slope close to and spectral levels 
decreasing by a factor 3 from the shoals to the middle shelf (dotted lines in figure 48b). 
The corresponding two-dimensional wavenumber spectra are not isotropic and show larger 
variance along a north-west to south-east axis, (e.g. figure 48a). Eor wavelengths under 300 
m, the topography is only resolved in the multibeam surveys conducted in small regions, in 
20 and 25 m, giving the two spectra and F^, respectively (figure 26), and represented 
as direction-integrated spectra in figure 48 (dashed lines). Other regions surveyed only 
along individual tracks, with the EM 1000 multibeam system or a ISIS 100 sonar (data ac¬ 
quired and provided by J. McNinch, ERE) gave unidirectional spectra comparable to these 
direction-integrated spectra. All spectra decrease approximately as (not shown), with 
variance levels between those of F® and F^■ 

In the model the topographic features that contribute to wave-bottom Bragg scat¬ 
tering (with horizontal scales smaller than five times the wavelength, see §III.B.4) are 
assumed to be statistically uniform on the entire shelf. A representative composite bottom 
elevation spectrum at these scales was generated from spectra computed from the 6” reso- 
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Figiue 47. Distribution of median grain sizes. 

Z) 5 o as a fimction of distance from the coast, for siu'frcial sediments in vibracore samples 
collected in the period 1994-1997 (Rebecca Beavers, Duke University, personal conmimii- 
cation, 1999) and Shipek grab samples collected dming SHOWEX (frgiue 35a). The dotted 
line indicates Dso = 0.15 mm, used in II. The average grain size Dso = 0.22 nun is indi¬ 
cated by the dash-dotted line, and the solid line is a polynomial frt used for hindcasts in the 
present chapter. 
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Figiue 48. Bottom elevation spectra 

(a) Composite two-dimensional bottom elevation spectnmi F^, based on 10-30 m depth 
bathymetiy in the entire model domain, and high-resolution multibeam bath 3 mietry in a 5 
km X 5 km region in 20 m depth. Contoin values are logjo with F^ in m'^rad”^, 

at 0.5 intervals. Circles indicate the bottom components that interact in 20 m depth with 
waves from the east with freqirencies 0.05 (irmer cucle), 0.12 (middle cucle) and 0.25 Hz 
(oirter cucle). Axes imits are reciprocal wavelengths 4/ (27t) and ly/(2 k). (b) Conespond- 
ing duectiou-integrated spectnrm (solid) and other spectra estimated from high-resohrtion 
mirltibeam bathymetry (dashed), and well resolved regions in the NOS database (dotted). 
The vertical lines mdicate the bottom scales responsible for scattering 0.08 Hz swell, for 
different incidence angles 0/. 
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lution bathymetry grid, ensemble averaging over regions with water depths between 10 and 
30 m (figure 48a, and solid line in figure 48b). The variance at large scales is lower than for 
the two spectra determined from regions well-resolved in the NOS database (dotted lines in 
figure 48b), possibly due to the coarse resolution of the bathymetry data in some regions. 
At higher wavenumbers, not resolved by the grid, the variance densities were taken from 
the small-scale spectrum with the larger variance in figure 26). Computations of the 
Bragg scattering source term using this spectrum may therefore overestimate scattering by 
the smaller scales (//27t > 0.002) and possibly underestimate scattering by the larger scales 
(l/ln < 0.002). To investigate this effect the model was also run with a similar composite 
spectrum in which the small scale part was replaced by (figure 26). 


5. Source terms 


The wave-bottom Bragg scattering source term Suragg is estimated using a uniform 
bottom elevation spectrum and integrated with a semi-implicit scheme, taking advantage of 
the linear relation between £’(/,.) and ^Bragg (/, ■) (§ III.C). Implementation of the bottom 
friction source term Smc is made more complex by the empirical parameterizations of the 
bottom Nikuradse roughness kj^. In the present chapter we generalize Tolman’s (1994) 
decomposition of in a ripple roughness kr and sheet flow roughness ks. 


kr 

ks 


Qb X Ai 



Ai 

1 


,2.8 


-0.4 


0.57 


[g(.-l)]i-" (271)2’ 


(V.2) 

(V.3) 


where and Ub are the r.m.s. amplitude of the bottom wave orbital displacement from the 
equilibrium position (half of the orbital diameter) and velocity at the top of the boundary 
layer, s is the sediment specific density and g is the gravity acceleration. While Tolman 
(1994) used values of the empirical coefficients Ai and A 2 determined by Madsen et al. 
(1990), we will examine if adjustments to these coefficients improve hindcasts of waves on 
the North Carolina shelf (figure 49). 

When the wave motion is not strong enough to generate vortex ripples, i.e. for val¬ 
ues of the Shields number less than a threshold xjTn-, kj\f is given by a relic ripple roughness 
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Te as a function of the Shields number v|/ noimalized by its critical value V|/c for sedi¬ 
ment motion. The original parameterization proposed by Tolman (1994) using Madsen 
et < 7 /. (1990) ripple roughness, is compared to the JONSWAP linear damping somce term, 
and an improved roughness parameterization is given by = 0.4, Ai — —2.5, ^3 = 1.2 
and ^4 = 0.05, and it is used below in the somce term Sfncj, giving better hindcast results, 
for the Shields number covered by the present data set (at X3: 0.03 < < 6.5). 
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ky:^. Madsen et al. (1990) proposed 


\\rry=A3\^c, 


(V.4) 


with A 3 = 1.2 determined empirically, and '\\fc the critical Shields number for the initia¬ 
tion of sediment motion on a flat bed and under sinusoidal waves. \|/c should be a known 
property of the sediment. In the absence of any data, the relic ripple roughness kyr was 
first proposed to be equal to the grain size diameter by Graber & Madsen (1988), which 
is appropriate for a flat bed, but Tolman (1994) suggested a larger value to account for 
relic bedforms and bioturbation, and proposed a constant value krr = 0.01 m. Our ear¬ 
lier hindcast study showed that this constant value yields good wave height predictions for 
very small waves (§ II.D.l), but the dissipation is so weak that different roughness values 
could yield similar wave heights. For larger forcing conditions (but still in the relic ripple 
regime), recent experiments suggest that kyy actually decreases with tjt on a smooth cohesive 
bed (Alex Babanin, University of South Australia, analysis of Lake George experimental 
data, personal communication, 2001). However hindcasts of waves on the North Carolina 
shelf (dominated by non-cohesive sand) presented in this chapter, suggest that the rough¬ 
ness may increase with the forcing, possibly due to the presence of relic bedforms that are 
more likely to be intercepted by larger orbital diameters. We therefore propose to use 

kyy = max { 0 . 01 m,A 4 a^} for \\f < \|/rr. (V.5) 


Al, A 2 , A 3 and A 4 define the response of the bed roughness to the wave forcing and the 
resulting wave energy dissipation, A 1 = 1.5, A 2 = —2.5, A 3 = 1.2 and A 4 = 0 corresponds 
to the original Tolman (1994) parameterization. 

In this parameterization, the critical value \\fc must be known beforehand. It is 
usually determined experimentally as a function of the sediment physical properties. Al¬ 
though recent experimental results for mixed sands suggest values of \|/c larger by about 
30% (Wallbridge et al., 1999), we use an analytical fit to Shields’ (1936) laboratory data: 


Vc 




0.3 


1 + 1.2D* 


+ 0.055 [1-exp (-0.02D*)], 


D 50 
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where V is the kinematic viscosity of water. This expression was proposed by Soulsby 
(1997), and includes the correction for very fine grain sizes by Soulsby and Whitehouse 
(1997) intended to limit \|/c to values less than 0.3. It is consistent with the laboratory data 
reviewed by Madsen & Grant (1976), and the ripple roughness parameterization of Madsen 
et al. (1990). 

The critical Shields number \|/c actually depends on the flow conditions and bed 
geometry. Knowledge of the exact threshold under irregular waves is not critical for the 
present work, and t|/c is better interpreted as a well-defined sediment property defined for a 
flat bed and given flow conditions (e.g. steady flow). In contrast, the parameterization of the 
equivalent ripple roughness kr, that also depends on the flow conditions and bed geometry, 
should be valid in field conditions. The laboratory experiments of Madsen et al. (1990) do 
not include all the properties of natural boundary layers, thereby giving a parameterization 
that may not be general enough. In particular it has been observed that ripples in the field 
are usually less sharp than ripples in the laboratory, perhaps as a result of the directional 
spreading of waves. This effect can be expected to reduce kr and may be parameterized by 
decreasing A i or increasing A 2 . The results of other local effects are not so clear, such as the 
impact of biological benthic activity (Fries et al., 1999). The detailed composition of mixed 
sands may also have some effect on the empirical parameters A 1,^2 and A 3 , and using the 
median grain diameter D^q in (V. 6 ) may not give the best value kyy for \\f ~ \[t].r should 
also be a function of the time varying properties of turbulence in the bottom boundary layer, 
as only the largest waves in a wave group may be able to generate vortices (responsible for 
the form drag and the very development of ripples, see Ayrton, 1910), but it should also 
depend on the past history of the bed morphology as relic ripples may determine if the 
vortices can be generated. 

In addition to these complications, due to the local properties of natural environ¬ 
ments, the practical use of bedform roughness parameterizations is hindered by spatial 
variations of the water depth and sediment grain size (see for example Green, 1986; fig¬ 
ure 35a), on scales not resolved by operational wave prediction models. The adaptation 
of the local parameterization to a larger scale (of the order of 0.1 to several kilometers). 
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was envisaged by Tolman (1995), and performed in chapter II using the known statistics 
of depth variations, and ignoring the variations of D 50 , because of insufficient data. This 
representation of subgrid variations in the properties of the bottom boundary layer is used 
here. Additional representation of subgrid variations of median grain sizes, assuming a 
Gaussian distribution, did not give significantly different results. 

A widely used alternative to this physics-based bottom friction source term assumes 
Sfric to be a linear function of the bottom velocity spectrum, with a proportionality coeffi¬ 
cient r/g^. r is an empirical constant with dimensions of m^s"^. This formulation emerged 
after the wave data analysis in the JOint North Sea WAve Project (JONSWAP), as varia¬ 
tions in the cross-shore energy flux were examined in terms of F, an attenuation coefficient. 
This linear dissipation can only be justified by assuming that the wave orbital velocities are 
weak compared to the mean (tidal) currents, but this should cause a modulation of F by 
the current that was not observed. However observations showed large variations of F 
(from 0.0019 to 0.160 m^s^^). Nevertheless this parameterization, with the averaged value 
F = 0.038 m^s"^ observed during JONSWAP, has encountered some success (Bouws and 
Komen, 1983), and replaced quadratic drag formulations proposed previously (Hasselmann 
and Collins, 1968). This ‘JONSWAP’ bottom friction source term, SfUc,;, may yield good 
results over sand because the corresponding dissipation factors fe decreases as a function 
of the Shields number, following the movable-bed model for relic roughness and sheet 
flow conditions (figure 49). However values of fe in the active ripple regime (\|/ > t|/rr) are 
significantly smaller than those observed over ripples, and Smcj fails to reproduce the ob¬ 
served strong attenuation of waves when active ripples are present (§ II.D.l). Fortunately 
for operational wave models, these active ripple conditions occur only in a small range of 
wave conditions, although these may be dominant in some regions. 

It must be emphasized that no parameterization has been tested for very large 
Shields numbers tl//tl/c > 5, for which very little laboratory or field data is available, but 
these conditions occur in large storms that are most important for navigation safety or 
coastal engineering. It is expected that a physics-based parameterization performs better 
than the JONSWAP empirical formulation calibrated in very different conditions. 
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5fric parameterizations used in the present chapter include 5fric,T, given by Tolman 
(1994) (i.e. Ai = 1.5, A 2 = —2.5, A 3 = 1.2 and A 4 = 0), and an improved movable bed 
parameterization Smcj for which the values of the empirical parameters have been adjusted 
to increase the overall hindcast skills: Ai = 0.4 and A 2 = —2.5, A 3 = 1.2, and A 4 = 0.05. 
Although no physical interpretation can be given for the JONSWAP source term Sfnc,;, it is 
used here for reference. 

D. MODEL-DATA COMPARISONS 

To objectively assess the accuracy with which the energy balance equation (V.l) 
describes the evolution of swell across the North Carolina continental shelf, a statistical 
analysis of hindcast results for the entire data set is presented here, based on 528 1-hour 
swell-dominated records (22 days) of SHOWEX observations and 121 three-hour swell- 
dominated records (15 days) of DUCK94 observations, out of the 2100 SHOWEX ob¬ 
servation records (87 days), and 725 DUCK94 records (91 days). Eewer records were 
available for some instruments, in particular A, B and WR(ERE) which failed during the 
1994 Hurricane Gordon, and X 6 which malfunctioned at the peak of the 1999 Hurricane 
Eloyd. The period in September and October 1999 when X 6 was not deployed was entirely 
excluded. 

The directional properties of the waves, influenced primarily by refraction and 
Bragg scattering (chapter III), are discussed first using only SHOWEX data. This is fol¬ 
lowed by an analysis of the wave heights, predominantly influenced by refraction and bot¬ 
tom friction using both DUCK94 and SHOWEX hindcasts. Eor all parameters the model 
errors generally grow across the shelf towards the shore. The ‘initial’ model errors, at the 
instruments located close to the offshore boundary, are discussed separately in § V.D.4. 

1. Mean direction 

A mean direction 0p at the peak frequency fp was computed for each instrument, 
using an energy-weighted average over a finite bandwidth of 0.15/^ centered at fp. Eol- 
lowing standard conventions (e.g. Kuik etal., 1988), 0p is defined as the direction of the 
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moment vector(ai, b\), with 


f‘\.015fp r27Z 

(ai.p,hi,p) = / / (cos0,sin0)£(/,0)d0d/ (V.8) 

J0.925fp Jo 

fp was determined from the measured frequency spectra at XI, so that the modeled and 
observed directions correspond to the exact same frequency band. Bp is well predicted 
by the model with no source terms or bottom friction only, with a maximum root mean 
square (r.m.s.) error of 8-10 degrees on the inner shelf, that decreases onshore as refrac¬ 
tion narrows the incoming spectra toward 70°, the beach-normal direction (figure 50). The 
conditions summarized here contain many events with large oblique angles that are well 
reproduced by the model, such as swell from Hurricane Floyd with offshore directions of 
160° refracted to 90° at 8M. The relatively large bias at 8M is not observed in similar hind- 
casts of DUCK94 swells (-1.6° bias, and r.m.s. error 6.9°) and may be related to a change 
in bottom topography not represented in the bathymetry grid used for ray computations. 
Errors at X5 are anomalously larger. In particular Bp at X5 was observed to oscillate be¬ 
tween 90 and 140°, whenever the offshore direction at X6 was slowly varying between 110 
and 130°, for all frequencies between 0.05 and 0.14 Hz. This suggests that the wave rays 
are bent by refraction, over bottom features not resolved in the sparse bathymetric surveys 
conducted on the outer shelf (Herbers et al, 2000). However, the errors are generally small 
compared to the changes in Bp across the shelf (up to 70° for the southerly swell from 
Hurricane Floyd), and are of the order of the model directional resolution (5°). 

The addition of Bragg scattering in the model biases the mean direction further to 
the north. This additional bias is due to the anisotropy of the bottom elevation spectrum 
used in the hindcasts, at wavelengths {In/I) less than 500 m, in combination with refraction 
by the large scale topography. These scattering bottom undulations have crests preferen¬ 
tially aligned with a south-west to north-east direction and have a strongest effect on waves 
propagating along that direction. The model bias and r.m.s. error added by Bragg scatter¬ 
ing is partially cancelled by bottom friction (compare runs 3 and 4 in figure 50). Indeed 
the scattered waves propagating at larger angles relative to the depth contours propagate on 
the shelf for a longer time, over which the cumulated effect of bottom friction is larger. A 
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8M WR(FRF) XI 


Bragg 

5. different bottom spectrum 

6. X6 only for b.c. 

X5 44014 X6 



Figiue 50. Model enors on the mean direction 0^ at the peak frequency 
(a) Model bias and (b) r.m.s. error for swell-dominated periods obser'ved dining SHOWEX. 
Rims 1^ use different sets of soince terms in (V.l). Soiuce terms in nms 5 and 6 are 
identical to those in nm 4, but nm 5 uses a different bottom elevation spectnmi that has 
a lower variance at small scales figiue 26b). Rim 6 is forced with X6 data at the 
offshore boimdary instead of the mterpolation of X6 and 44014 data used in all other nms. 
Positive values in (a) correspond to a clockwise bias. 
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weaker variance in the bottom elevation spectrum for In/l < 500 m clearly reduces this 
effect (run 5 in figure 50). 

2. Directional spread 

The directional spread at the peak frequency Oe ^ was computed for each instru¬ 
ment, as was the mean direction Op, for a finite frequency band around the peak of the 
measured frequency spectrum at XL Oe^p, defined in radians, as 

/■1-075/p flTt 

Oe,p = [2(l-ai,pCos0p-Z7i,pSin0p)]i/V / / £(/,0)d0d/, (V.9) 

./ 0 . 925 /p Jo 

is an estimate of the standard deviation of the directional distribution of wave energy, 
and thus can be loosely interpreted as the half-width of the directional spectrum. For an 
isotropic spectrum Oe p is 2^/^ radians, that is 81°. 

Observed directional spreads are, in general, relatively uniform across the shelf, de¬ 
creasing slightly towards the coast. Observations at X6 range from 10 to 55° with typical 
values between 30 and 40°, whereas observations at 8M vary between 12 and 27° with typ¬ 
ical values of 15-25° (figure 51). However, when the offshore directional spectrum is very 
narrow (oe,p (X6) < 20°), Oe^p increases towards the shore. This observation is contrary 
to the general belief that waves become one-directional in shallow water, as expected from 
depth-refraction. Indeed model calculations that incorporate refraction, as well as bottom 
friction (in order to provide reasonable wave heights) give directional spreads much nar¬ 
rower than the observations (figure 51). This bias, largest at 8M and on the inner shelf, can 
be observed throughout the shelf. 

Statistics from model runs with different source terms (figure 52) demonstrate how 
bottom friction and wave-bottom Bragg scattering contribute to Oe p. Run 1, without source 
terms, clearly underestimates Oe p and only contains the narrowing effect of refraction 
as waves approach the shore. Run 2 includes bottom friction that narrows the spectrum 
further: waves propagating at larger angles relative to the depth contours are more exposed 
to dissipation, because of their longer propagation time across the shelf. Wave-bottom 
Bragg scattering opposes these two narrowing effects by diffusing the energy around the 
mean direction (see § III.C). While Bragg scattering alone tends to give spectra that are too 
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Measured ao p at X6 

Figure 51. Nearshore versus offshore directional spread 
Diiectional spread at the peak frequency at 8M, as a fimction of the offshore diiectional 
spread measiued at X6. Tlie solid Ime separates spectia that are broader in the nearshore 
and spectra that are broader offshore. 
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Figiue 52. Model eiiors on the diiectional spread ae,^ at the peak frequency 
Same foimat as frgiue 50. The r.m.s. eiior is replaced by a scatter index defrned as the ratio 
of the r.m.s. eiior and the r.m.s. value 
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broad in the middle of the shelf (figure 52, run 3), the addition of Bragg seattering to bottom 
frietion generally yields good agreement aeross the shelf (figure 52, run 4), redueing the a 
seatter index from 0.51 (run 2) to 0.25 (run 4) at 8M, with eomparable improvement on the 
inner shelf (figure 52). 

However, results are sensitive to the ehoiee of the bottom elevation speetrum. Run 
5 using a different speetrum with the small-seale varianee taken from F 2 ^ (figure 26b, 
representative of smoother regions of the shelf) only partially aeeounts for the observed 
inerease of Oe,p (figure 52, run 5). Bragg seattering most strongly affeets narrow speetra, 
and reproduees well the observed broadening of very narrow offshore speetra (figure 51). 
For broad offshore speetra (20 < OQ^p (X6) < 40) the observed speetra at 8M are still 5- 
10° broader than the predieted speetra. The predieted seattering effeet is also strongest for 
waves with low peak frequeneies (figure 53), but the narrowest offshore speetra tend to 
eoineide with the ones with the lowest peak frequeney, so that the influenee of the speetral 
width and peak frequeney eannot be elearly separated. 

As the effeet of Bragg seattering inereases with the bottom elevation varianee, it is 
likely that the bottom speetrum used in run 4, with relatively large small-seale varianee, 
represents an upper bound on the true effeet of wave-bottom Bragg seattering aeross this 
eontinental shelf. The remaining negative bias of Oe p predietions in run 4 suggests that 
other proeesses are important. Higher order Bragg seattering, among others, may be the 
eause of further direetional spreading. 

The bottom frietion souree term was somewhat arbitrarily assumed to be isotropie. 
A quadratie drag law gives a differenee between the major and minor prineipal axes of 
the direetional dissipation souree term that is up to a faetor 2 for unidireetional waves 
(Hasselmann & Collins, 1968). This type of anisotropy would eause more dissipation for 
waves with direetions along the mean direetion, and less dissipation in the perpendieular 
direetions, therefore inereasing the direetional spread. The absenee of wave refleetion from 
the beaeh in the model is another likely eause for the negative bias. Even though the 
offshore traveling waves are removed in the direetional speetrum estimates and Oe p at 8M, 
they eontribute to the Oe p estimates at other sensors. 
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Figure 53. Relative error on the directional spread Oe^p as a function of the peak frequency 
Relative errors of 1 or 0.5 represent an overprediction or underprediction by a factor 2, 
respectively. 
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3. Wave heights 

For the analysis of wave height predietions the SHOWEX data set is augmented 
with the DUCK94 data. Statisties given below eombine results at SHOWEX and DUCK94 
instruments loeated nearby or in similar water depths. Changes in the peak frequeney 
aeross the shelf are negligible, but the signifieant wave height Hs generally deereases from 
offshore to the nearshore with observed Hs at 8M about half of the observed values at 
X6 for moderate-energetie swell, and a generally smaller reduetion 25%) in benign 
eonditions (Hs (X6) < 1.5m). This effeet is explained only in part by refraetion that reduees 
the wave heights of waves that propagate onshore at large oblique angles relative to the 
depth eontours. Indeed the model without souree terms (figure 54, run 1), that aeeounts 
only for the effeets of refraetion and shoaling, overprediets wave heights with a typieal bias 
of 0.2 m on the inner shelf, and gives an overall (for all sensors) seatter index of 0.26 for Hs. 
Adding Bragg seattering slightly inereases model errors: the inereased direetional spread 
translates into a larger eross-shelf propagation time (on average), whieh, in the absenee of 
dissipation, inereases the wave height, giving an overall seatter index of 0.29 (figure 54, 
run 3). 

Ineluding bottom frietion dramatieally reduees the model errors. Run 4b, using the 
Tolman movable-bed souree term, based on laboratory data without any empirieal tuning 
to field eonditions, reduees the overall seatter index of Hs predietions from 0.29 to only 
0.15 (figure 54). This result supports the hypothesis that the formation of vortex ripples 
(observed in ehapter IV) and their feedbaek on the waves through enhaneed bottom rough¬ 
ness, are the primary meehanisms for for wave attenuation aeross a sandy eontinental shelf. 
However, a negative bias of -9 em at 8M indieates a model tendeney to overestimate dissi¬ 
pation, already noted in ehapter II. 

On average the empirieal JONSWAP bottom frietion souree term performs about 
equally well (figure 54, run 4a), giving an overall seatter index of 0.16. However, this 
souree term gives poor results at CHEV2, a site loeated down wave of shallow shoals (but 
not shallow enough for depth-indueed breaking), with a seatter index of 0.53 and a positive 
bias of 20 em. This bias is the result of large model-data diserepaneies during the arrival 
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Figme 54. Model eiTors in predictions of the significant wave heights 
(a) Model bias and (b) scatter index (ratio of r.m.s. enor and r.m.s. value) for H^, m 
swell-dominated periods obseiwed dining SHOWEX and DUCK94. Rim 1 includes no 
soiuce terms (refiaction and shoaling only), iim 3 includes wave-bottom Bragg scattering, 
and different friction soiuce teims are added in nms 4a (JONSWAP) 4b (Tolman) and 4c 
(improved Tohnan). Soiuce terms in nm 6 are identical to those in nm 4c, but only X6 
data is used for the offshore boimdary condition in SHOWEX cases instead of a linear 
interpolation of X6 and 44014 data. For DUCK94 cases the offshore boimdary condition 
is obtained fiorn 44014. 
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of swell from Hurricane Gert (17-21 September 1999). Observed wave heights at CHLV2 
during this event are 70% smaller than predicted by refraction and shoaling, which corre¬ 
sponds to a dissipation of 84% of the incident wave energy flux. This strong dissipation 
is probably the result of active ripple generation on the shoals offshore of CHLV2 (figure 
42), and the JONS WAP source term grossly underpredicts bottom friction in active ripple 
conditions (figure 49). The improved Tolman source term Smcj (using a modification of 
Madsen et a/.’s (1990) empirical parameters, designed to reduce the model errors) gives an 
overall scatter index of 0.13, with a maximum of 0.20 at X2 (run 4c in figure 54). 

The distribution of model errors as a function of the Shields number, estimated 
from waves measured at representative inner shelf sites X3 (SHOWEX) and C (DUCK94), 
is shown in figure 55. The largest errors of a non-dissipative model (run 1) clearly occur 
for large values of the Shields number when the presence of active ripples is expected. A 
model using the JONSWAP friction source term also gives the largest errors in this regime. 
It should be noted that the two largest Shields numbers correspond to DUCK94 3-hour 
records in the afternoon of 18 November 1994, after the peak of Hurricane Gordon (see 
figure 16), as the eye hurricane was moving back to the south (figure 43). This is the only 
data available for large hurricane-swell forcing conditions, with near-sheet flow conditions 
(the sheet flow conditions indicated in figure 18 correspond to the peak of the storm a few 
hours earlier, on the morning of 18 November). In the strongest forcing conditions the 
JONSWAP source term largely overestimates the wave height. As the dissipation factor 
is expected to increase with the onset of sheet flow, the JONSWAP source term, with a 
decreasing dissipation factor, may systematically overpredict wave heights in sheet flow 
conditions. 

In contrast the standard Tolman source term Smcj significantly underestimates 
wave height in active ripple conditions (figure 55g-i). If the ripples (observed on most 
of the shelf, chapter IV) had been as rough than those in Madsen’s (1990) laboratory ex¬ 
periments, 5fric,T should give reasonable wave heights. The observed discrepancies suggest 
that ripples in the field are significantly smoother than in the laboratory. This may be due 
to the directional spreading of natural waves, as well as the initial bed configuration. Both 
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Figme 55. Relative eiior ofHs predictions as a fimctiou of the Shields uiunber. 
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may increase the shape of the ripple crests and the number of defects in the ripple pattern, 
and modify the equilibrium bed configuration. A reduction of the maximum ripple rough¬ 
ness by a factor 4 gives better agreement with the DUCK94 and SHOWEX data (run 4c, 
using 5 mc,i, figure 55j-l). 

The present dataset is too limited to address effects of past history of the bedforms, 
and there is no clear indication that the ripple roughness should be different between wave 
events with increasing, stable, or decreasing forcing conditions. The behavior of Smcj and 
5fric,i seems acceptable in the strong forcing conditions prevailing during the Hurricane 
Gordon event, although more data would be needed for a more quantitative assessment. 

4. Model errors and offshore boundary 

Errors are not zero at the buoy (X6 and/or 44014) used to provide the offshore 
boundary condition, as only incoming waves are used for the forcing whereas the outgoing 
waves are radiated by the model. Errors are also introduced by the numerical propagation 
and interpolation between the buoy and the grid boundary. 

Errors at X6, when the boundary condition is determined from X6 data only (run 
6 ), are entirely the result of the difference between outgoing waves predicted by the model 
and caused by refraction and back-scattering, and those estimated from the buoy data us¬ 
ing the Maximum Entropy Method. Since this method is only constrained to fit the first 
four Eourier components of the directional spectra, it may give significant errors for these 
outgoing waves that usually carry very little energy, compared to the incoming waves. The 
corresponding r.m.s. error on the mean direction is 5°, and the scatter indices for Oe ^ and 
Hs are 0.18 and 0.06, respectively. 

Eor this same run 6, errors at 44014 are enhanced (11° r.m.s. for Op, and scatter 
indices for Oe.p and of 0.33 and 0.17, respectively) possibly due to propagation errors 
across the shelf break, but also likely caused by spatial variations in the offshore wave field, 
over the 60 km separating X6 and 44014. These offshore variations may be caused in part 
by wave-current refraction in the Gulf Stream and its eddies (Holthuijsen & Tolman, 1991), 
often seen along the shelf break in sea surface temperature satellite images. This effect 
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should be strongest at high frequencies and could explain the fact that the energy difference 
between X6 and 44014 for remotely generated swells is either constant or increases with 
increasing frequency. The errors on the directional spread Oe p may also be influenced by 
differences in the response between a Waverider buoy (X6) and a 3-m discus pitch and roll 
buoy (44014), the latter giving generally larger directional spreads (O’Reilly et al, 1996). 

These errors in run 6 give an order of magnitude for the model errors due to the 
treatment of the boundary condition in all other runs. The imperfect knowledge of the 
variations in the wave field along the shelf break is probably a major source of errors in 
runs 4a-c. Based on errors at 44014 and X6, contribution to the scatter indices resulting 
from boundary condition errors is estimated to about 0.1 for Hs and 0.15 for Oe,p. In view 
of this high background error in the treatment of boundary conditions, the energy balance 
(V.l) gives a good description of swell evolution. 

E. SUMMARY AND PRACTICAL IMPLICATIONS 

The numerical model CREST was implemented on a large portion of the North 
Carolina-Virginia continental shelf for a comprehensive hindcast of all swell-dominated 
conditions observed during the experiments DUCK94 and SHOWEX. The model inte¬ 
grates the energy balance (V.l) from deep water (1500 m depth) to the inner shelf (8 m 
depth), with little numerical diffusion, allowing detailed comparison of measured wave 
parameters with results obtained using various sets of source terms. The effect of the 
small-scale (comparable to the surface wavelength) shelf topography, was represented by 
a wave-bottom Bragg scattering source term, assuming uniform statistics determined from 
high-resolution bathymetric surveys, explains most of the observed broadening of the wave 
spectrum towards the shore, that occasionally balances the narrowing caused by refraction 
over the quasi-plane large-scale bathymetry. Predicted directional spectra are nevertheless 
still too narrow on the inner shelf (by about 20% compared to 50% for hindcasts without 
Bragg scattering), maybe as a result of wave reflection from the beach, anisotropy in the 
bottom friction source term, and/or higher order wave-bottom Bragg scattering, that are not 
represented in the present model. The strong attenuation of large swells (inferred dissipa- 
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tion up to 84% of the incident wave energy flux) is well reproduced by Tolman’s (1994) 
bottom friction source term, that accounts for the generation of sand ripples by waves and 
their feedback on the waves. Including this source term into the spectral refraction model 
reduces the scatter index for the significant wave heights Hs from 0.27 to 0.15. A large 
fraction of the remaining error may be due to the imperfect knowledge of the offshore 
boundary condition which may cause a scatter of 0.1 on the outer shelf, in predictions of 
Hs by a perfect model. 

The widely used JONSWAP source term gives slightly larger errors on average 
(overall scatter index equal to 0.16), but occasionally causes much larger errors than the 
Tolman source term, because it does not account for large increases in bottom roughness 
under active ripple conditions. Tolman’s (1994) movable bed model tends to overestimate 
dissipation in strong forcing conditions, suggesting that ripples in the field are significantly 
smoother than in Madsen et a/.’s (1990) laboratory experiments, that were used to calibrate 
the empirical ripple roughness. The reason for this difference between laboratory and field 
conditions is not clear but may be caused by the universal presence of relic bedforms in 
the field, and the directional spreading of natural waves. An ad hoc reduction of the ripple 
roughness by a factor 4 improves the hindcast results and reduces the overall scatter index 
for the wave heights from 0.15, for Tolman’s (1994) source term, to 0.13. 

Overall the energy balance (V.l) with a bottom friction source term, dissipating 
swell energy, and a wave-bottom Bragg scattering source term (chapter III), broadening 
the directional spectrum in the cases presented here, provides a good description of swell 
evolution across the North Carolina continental shelf. In addition to the direct local effects 
of each source term, the numerical integrations of (V.l) reveal indirect effects coupling 
the source terms as waves propagate over a finite distance. Bottom friction reduces the 
directional spread since obliquely traveling waves are more exposed to dissipation as they 
propagate for a longer distance (and time) across the shelf, while Bragg scattering lengthens 
the average propagation distance. In the absence of dissipation this longer propagation may 
yield larger wave heights, in a way similar to shoaling waves propagating straight towards 
the shore: the onshore energy flux is conserved, but the average onshore group speed is 
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reduced by scattering into obliquely traveling components, so that the wave height must 
increase. However, in the presence of dissipation the opposite happens: the wave height is 
reduced because these obliquely propagating waves have more time to be dissipated. 
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VI 


CONCLUSIONS AND PERSPECTIVES 


The work presented in the previous chapters focused on the physical processes and 
numerical prediction of the evolution of swell in shallow water, using new high-quality 
field data from the North Carolina continental shelf. 

A new numerical wave model named CREST (Coupled Rays with Eulerian Source 
Terms) was presented in chapter II, that integrates in time an energy balance equation for 
the evolution of the wave frequency-directional spectrum. The novelty of this model resides 
in its hybrid numerical scheme. It couples inverse ray trajectories adopted from Eagrangian 
models, tracing back the paths followed by wave groups, and a modification of the wave 
energy along these paths using a spectral source term. The source term, evaluated on the 
Eulerian grid from which the rays are traced, and interpolated onto these rays, prescribes 
how the energy of each component changes in time based on the entire wave spectrum 
and environmental factors such as the bottom sediment nature and the wind speed. This 
hybrid method essentially eliminates the numerical diffusion of finite difference schemes, 
allowing accurate and high-resolution predictions of the wave spectrum. We demonstrated 
that this numerical method can be applied to a large area (100 km x 400 km for the North 
Carolina-Virginia continental shelf) and implemented it on a modern workstation, because 
the grid, on which source terms are evaluated, can be much coarser than the bathymetry 
grid used for the ray computations. Although detailed comparisons have still to be made, 
it can be expected that the computational cost of the Eagrangian advection scheme in its 
present form is comparable to or larger than the cost of a finite-difference scheme, which 
uses higher resolution grids to resolve scales responsible for depth refraction. However 
the new model is expected to be more efficient for the prediction of the generation and 
non-linear evolution of wind sea spectra, because the costly source term computation in 
finite difference models is performed on the same high-resolution grid as the advection, 
while it is done on a coarse source term grid in the present model. CREST may therefore 
be used as a research community tool for hindcasting wind seas with exact calculations of 
the quartet wave-wave interaction source term. The computational cost of CREST may be 
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reduced by adopting so-called ‘piecewise ray’ methods that use rays ray segments over sin¬ 
gle time steps (e.g. the TOMAWAC model, Benoit et al, 1996). Practical trade-off criteria 
should be determined for choosing intermediate schemes with ray advection over several 
time steps, reducing the numerical cost at the expense of numerical diffusion. These cost¬ 
saving measures may promote an operational or engineering use of the code as computers 
keep getting faster. However, the hybrid Eulerian-Lagragian scheme, in its present form, 
cannot handle time-varying depth and strong variable currents that are important for wave 
evolution in macrotidal coastal regions (e.g. the North Sea, the English Channel, or the 
Erench Atlantic continental shelf), where ray trajectories vary in time. 

Complementing earlier wave observations made on the North Carolina continental 
shelf during the DUCK94 experiment, a new experiment, SHOWEX, provided high-quality 
directional wave measurements across the same region, over a three-month period with a 
wide range of wave conditions. CREST was implemented on the North Carolina conti¬ 
nental shelf and produced hindcasts of swell observed during both experiments, used to 
quantitatively evaluate wave-bottom interactions in the absence of currents. The evolution 
of swell in moderately shallow water is essentially influenced by the bottom topography at 
all scales. The field data and model results confirmed that large features with horizontal 
dimensions much larger than the wavelength cause well-known refraction and shoaling ef¬ 
fects, such as increased wave heights around headlands, and provided new insights into the 
effects of smaller scale bottom topography and roughness. 

Bottom undulations with horizontal scales comparable to the surface wavelenghts 
cause directional scattering of waves through a Bragg resonance. In a spectral description 
this wave-bottom interaction can be represented, at the lowest order in the bottom and sur¬ 
face slopes, as interactions among triads of two wave components with the same frequency, 
exchanging energy, and a bed undulation with the difference vector wavenumber that makes 
this exchange possible. The evolution of the wave spectrum caused by this triad interac¬ 
tion was derived in the form of a spectral source term Ssragg by Hasselmann (1966), using 
a Hamiltonian formalism, and Eong (1973) showed that it could cause significant back- 
scattering (i.e. reflection) of swell towards deeper waters for a realistic bottom elevation 
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spectrum. Here the theory was extended to slowly varying random waves and bottom ele¬ 
vations on a sloping shelf, and a correction of the coupling coefficient was given. Hindcasts 
of swell transformation across the North Carolina shelf, marked by gentle sand ridges with 
spacings of 1 km or more, showed that back-scattering is insignificant, but ‘forward scatter¬ 
ing’ (the interacting waves have almost the same direction) increases the directional spread 
of the wave spectrum. This unexpected effect, not considered in previous studies, op¬ 
poses directional narrowing resulting from depth refraction. The combination of these two 
processes in CREST, a numerical model with little numerical diffusion, yields directional 
spread hindcasts that agree well with directional wave measurements from SHOWEX. In 
8-15 m depth the use of the Bragg scattering source term increases the predicted spread by 
a factor 2-3, in particular for low frequency swell. Applications to other regions where de¬ 
tailed bathymetric information and directional wave data is available should provide further 
verification of the wave-bottom Bragg scattering theory for random waves and topography. 
Such studies should encourage new work on the practical implications of higher order the¬ 
ories and effects of currents on the wave-bottom interactions. Wave-bottom Bragg scatter¬ 
ing is the simplest resonant interaction involving surface gravity waves, but it shares some 
properties of the non-linear wave-wave interactions. Thus the ideas developed in chapter 
ITT on the detuning effects of the bottom slope on the resonant interactions may be applied 
to wave-wave interactions, providing practical guidelines for the use of resonant interaction 
theory in varying depth. 

At scales smaller than the wave orbital motion, sand ripples and biogenic mounds 
on the seabed contribute to the bed roughness experienced by the waves. It is well-known 
that non-cohesive sediments are mobilized by sufficiently energetic wave conditions form¬ 
ing ‘vortex ripples’ that dramatically enhance dissipation of the wave energy. Widespread 
ripple marks were observed across the North Carolina using sidescan sonar (chapter IV). 
The combination of sonar surveys and time-series of wave observations established their 
‘vortex ripple’ nature, consistent recent observations in shallower water. A parameteri¬ 
zation for the dissipation source term Sfnc of wave energy by bottom friction over these 
bedforms was proposed by Madsen et a/.(1990), based on careful laboratory experiments. 
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It was tested for the first time against field data from the DUCK94 experiment (chapter II). 
The model hindcasts reproduces well the observed swell attenuation across the shelf, but 
Madsen et a/.’s (1990) parameterization tends to overestimate dissipation in ripple-forming 
conditions. 

A comprehensive analysis of all swell-dominated conditions observed during 
DUCK94 and SHOWEX was conducted to objectively examine the performance of bot¬ 
tom friction parameterizations, and propose empirical adjustments (chapter V). The effect 
of Bragg scattering was included in the source term together with various bottom friction 
parameterizations. The widely used empirical ‘JONSWAP’ parameterization of the bot¬ 
tom friction gives average errors similar to Tolman’s (1994) source term that uses Madsen 
et aUs (1990) ripple roughness. However the JONSWAP source term grossly underesti¬ 
mates the swell decay in active ripple conditions. 

A modification of Tolman’s bottom friction source term, increasing the roughness 
of relic bedforms as a function of the wave forcing, and reducing the roughness of actively 
generated vortex ripples by a factor 4, improves the hindcast accuracy. These results sug¬ 
gest that natural ripples formed by a directionally spread wave field are less rough than 
those generated by uni-directional waves in wave flumes. Both Madsen etaUs (1990) 
parameterization and the one presented here are based on dissipation rates inferred from 
observed wave attenuation. A true verification of the parameterization requires detailed 
measurements of the properties of the bottom boundary layer under waves in realistic con¬ 
ditions. Reynolds stress profiles measured by T. Stanton & E. Thornton during SHOWEX, 
should help constrain parameterizations of the ripple roughness and the directional proper¬ 
ties of the bottom friction source term. 

In addition to the direct local effects of bottom friction (dissipating swell energy) 
and Bragg scattering (broadening the directional spectrum), indirect effects couple these 
source terms as waves propagate over a finite distance. Bottom friction reduces the di¬ 
rectional spread since obliquely travelling waves are more exposed to dissipation as they 
propagate for a longer distance (and time) across the shelf, while Bragg scattering lengthens 
the average propagation distance. In the absence of dissipation this longer propagation may 
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yield larger wave heights. However, in the presence of dissipation the opposite happens: 
the wave height is reduced because these obliquely propagating waves have more time to 
be dissipated. This interaction of the source terms in the energy balance equation shows 
that one source term cannot be considered in isolation, unless all others are negligible. It 
also emphasizes the importance of the directional distribution of wave energy for nearshore 
wave height predictions. 

Theories and parameterizations for wave-bottom interactions were examined here 
only for waves in the swell frequency band at a single site. The same interactions however 
also affect waves in the infragravity band (0.01-0.05 Hz) and the wind sea band (0.1- 
1 Hz). Extension and verification of the present work into these other frequency domains 
and other shallow water regions is needed for a comprehensive validation of Smc and Ssragg 
formulations. Incorporating these validated source terms in the wind sea energy balance 
may help improve parameterizations of other source terms: wind input, and dissipation due 
to deep-water wave breaking. In the absence of wind forcing, the present work establishes 
that the transformation of swell across a sandy continental shelf is well described by the 
energy balance equation 

= 5mc (k) + SBragg (k) , (VI. 1) 

where the bottom friction source term parameterizes movable bed processes: relic rough¬ 
ness, vortex ripple formation, and sheet flow. On the wide and shallow shelf of North Car¬ 
olina, these effects are locally weak: energy dissipation is of the order of 0.1-10 W m ^ 
compared to a typical energy flux of 10-100 kW m but they add up to a strong attenu¬ 
ation as waves propagate across the shelf, occasionally reducing the wave heights of large 
swells by more than 50 % . The formation of vortex ripples is therefore among the most 
important processes for the nearshore wave climate. 
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Figure 56. Singing sands, lake Huron, 28 April 2001 
Sand ripples, waves, and eapillary ‘Wilton ripples’. From auto-organization to resonant 
wave-wave interaetions. Photo by Fanny Girard-Ardhuin. 
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APPENDIX A. DERIVATION OF THE ENERGY 

BALANCE (III.47) 


1 . DERIVATION OF £2,2 (K) 

The governing equation (111.34) for 4>2(t) is an undamped forced harmonic os¬ 
cillator with a resonant frequency (O given by the dispersion relation (111.28). Applying 
a Fourier decomposition to the right hand side forcing terms, (III.34) can be written as a 
linear superposition of equations of the type 

^+(oVl=e‘“^ (A.l) 


In order to specify a unique solution to (A.l), initial conditions must be prescribed. In the 
limit of large propagation distances the initial conditions contribute a negligible bounded 
term to the solution. Following Hasselmann (1962), we chose f\ (0) = 0 and d/i/dt (0) = 
0, giving the solution 


/i ((0,(o';t) 


eico't _ gicot I (-(jj _ j (Q 

0)2 - (0'2 


for ^ d, 


(A.2) 


/i (to, mV) 

^>2 (t) is given by 


sin (tot) 
2m' 2i(o'(o 


for d 


±(D. 


(A.3) 


^2.k (0 = Ea (k,k') 5k_k'$i,k'/i {d -sd;t) , (A.4) 

k' 

where 

A(k,k') = [k —d^tanh{kH)] ^ ^ . (A.5) 

K 

The third-order energies £' 1.2 and £ 2,1 that result from the correlations between 
4>| ^k ^2 l^k (A-4) are found to be bounded but the fourth order energy £ 2,2 grows with 
time. Substituting (A.4) in (III. 16) and taking the limit (III. 17) to a continuous spectrum 
yields the solution for £ 2.2 (t,k) defined by (III.25): 

foo p2n ^ ^ 

£ 2 , 2 (Tk)=/ / A(k,k')£®(k-k')£2(t,k')|/i((O,-(o';t)rclk'd0' + ..., (A.6) 

./ 0 *-'0 
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where k' = k' (cos 0', sin0'), and the omitted terms (...) include bounded terms forced by 
surface non-linearity. Assuming that the frequency spectrum is continuous, the contribution 
of exact resonant interactions ik' = k) to (A.6) is negligible compared to contributions of 
near-resonant interactions {k' ~ k) that span a finite range of wavenumbers, and thus (A.2) 
can be substituted in (A.6). We now change the integration variable k' to (o', given by the 
dispersion relation (III.28), and replace |/i by the expansion in powers of ((O — (o'). 


^ ^ 2-2ccs|(».-».0>l+O(a,-».-) 

' ' '' 4M2((o-oy)- 

The leading term in this expansion is the real part of a function of the complex variable 
z = ((D —(o'), (l — e^^') / ( 2 ( 0 ^^), which is differetiable with respect to z on the entire 
Argan-Cauchy plane, except at z = 0. The integral /q°° |/i p dto' cuts through that singularity 
and its value nt/ (2(0^) can be obtained by contour integration on the complex plane. Other 
terms in (A.6) are continuous so that they can be considered constant in a small region 
around the singularity and we obtain 


^2,2 (Tk) = t 


2n 


47t(OA:'^ cos^ (0 — 0') 


sinh {IkH) [2kH -T sinh {2kH)] 
-l-bounded terms , 


(k-k')£2(Tk') d0' 


(A.8) 


For large t the derivative of £ 2,2 (T k) with respect to t yields equation (III.35). 


2. DERIVATION OF (K) + Ef^ (K) 

The particular solution to (III.37)-(III.39) in the vicinity of x = 0 can be written 
as 

(^f (x,z,0 = 1^ 

k,s 

The solution for the bound component 4>3follows from substituting the second-order 
velocity potential (III.32) in the bottom boundary condition (III.38 with term V set to zero) 

(0 = - E *^5k-k'^2,k (0 • (A. 10) 

k' 


^3,kW 


cosh {kz + kH) 
cosh {kH) 


+ 4 ) 3 ’® 


’kW 


sinh {kz + kH) 
cosh {kH) 


,ikx 


(A.9) 
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Substitution of (A. 10) in the surface boundary condition (111.39 with the right-hand side set 
to zero) yields the forced harmonic oscillator equation 


dt^ 




k k' 


-5, 


k k' 




(A.ll) 


Using (111.34) and (A.4) we have 


^ + (dM k = -Ea (k,k') 5k_k'E a (k',k") 5k^_k''^j,k"/i 

h-ff ’ 


k k' 


E tanh (kH) 5k_k' E^ (k'. k") 5k._k"^i,k" (0, 


(A.12) 


where A is defined by (A.5). The only terms of (A.12) that force growing correlations with 
(|)i are those that have a product of two factors A and equal wavenumbers k" = k. Other 
terms force bounded correlations and can be neglected. Therefore (A.12) can be regarded 
as a linear combination of forced oscillator equations of the form 

dVz 


dt^ 


+ (0^2 = /l ((o',-5(OT) • 


(A.13) 


The solution of (A.13) for ^ (Or and initial conditions /2 (0) = df 2 /dt (0) = 0 is 

. . ; X te‘®® - sin (rot) /ro (ro' -h ro) e^“'‘ -h (ro' - ro) - 2ro'e^“^ .. ,. x 

/ 2 (ro,ro',^;t)=—^ - , A -• (A.14) 


2isro(ro'2-ro2) 
Following the steps of § A. 1, we obtain 


2ro' (ro'2 - 0)2 


sc / SC / f^^An(i)k'^cos^(Q-Q')F^(k-k') 


Jo sinh {2kH) [2kH + sinh {2kH)] 
+ bounded terms. 


(A.15) 


Taking the derivative of (A.15) with respect to t yields equation (III.41). 


3. DERIVATION OF (K) + (K 

In the vicinity of x = 0, t])"® can be written in the form 


k,5 


cosh (kH) 
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where ^>3 satisfies 


(^ + ®^)^'ik = 2is(0^^e (A.16) 

With the solution given by eq. (A.3): 

<l>3,k(<)=-<-^e-“'. (A.17) 

This solution is correlated with the first order velocity potential, giving the energy contri¬ 
bution at X = 0 

dE2(Lk) 

(t, k) + 3 (t, k) = —t - —- + bounded term. (A. 18) 

Taking the derivative of (A. 18) with respect to t yields (111.42). 




and the remaining two terms follow from the bottom and surface boundary conditions. 
Substituting (A.20)-(A.23) in the bottom boundary condition (111.38, with III and IV set to 
zero) gives 

(A.24) 

Substituting (A.20)-(A.24) in the surface boundary condition (III.39 with VI set to zero) 
yields a forced harmonic oscillator equation for 4>3 (t): 


^ + ® )^'ikW=isge 
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+-k-Vkr + k-VH\^\^. 

K 


(A.25) 


(A.25) is of the same form as (A.l) with only resonant forcing terms ((£)' = ±(0) and a 
solution given by (A. 3). The covariance of and the first order potential, defined by 
(III. 16), is given by 

a 




2k 

+t^k- 
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^//tanh {kH) + + ktanh {kH) VH 
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t(0 
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k^ tanh {kH) 


where C„ is the group speed 


C =- 

^ kr 


1 


■ + 


krH 


(A.27) 


2 sinh(2k^//) 

and the bounded term is given by the initial conditions (the second right hand side term in 
A.3). From the dispersion relation (III.28) we find 


2k^ 


(A.28) 


and 


V[Cg tanh {KH)] _ CO 
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2k//tanh {kH) j- Vk^. 


(A.29) 



Using (A.28),(A.29), and the definition of the Lagrangian energy spectrum E 2 (kr) (111.45), 
(A.26) reduces to 




(k) 


Akr^^Akr^y = -tV ■ {CgE2 (k;. 




bounded terms (A.30) 


Writing Cg as (CgCOS0,CgSin0), where 0 is the local ray direction, the divergence term in 
(A.30) can be expressed in terms of intrinsic coordinates: 


V ■ [CgE2 (k^) Akr,x^kr,y\ = 


d [CgE2 (k^) Akr^x^r, 
dr 


+ CgE2 (k^) Akr^x^kr^y ^ 


dn 


(A.31) 


where r and n are the local along-ray and ray-normal coordinates. Longuet-Higgins (1957, 
equations 6, 10 and 21) showed that (A.31) can be simplified using ray theory. Defining 
the small separation p of two rays that are parallel in the vicinity of x = 0, we have 


do 1 dp 

dn p dr 


(A.32) 


and 


^ {pCgAkr^x^krj) 
dr 


= 0 


(A.33) 


Substituting (A.32) and (A.33) in (A.31), and taking the limit jak^l —^ 0, we have 


(k) -f £’i '®3 (k) = + bounded term . 


(A.34) 


Finally the derivative of (A.34) with respect to t yields (III.43). 
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